Milo on disease VS healthy dataset

suppressPackageStartupMessages({
  library(tidyverse)
  library(irlba)
  library(DropletUtils)
  library(scater)
  library(scran)
  library(Seurat) ## just 4 loading the object
  library(SingleCellExperiment)
  library(miloR)
  })

Using data from Ramachandran et al. 2019 (GEO accessiion: GSE136103). The Seurat object was downloaded from https://datashare.is.ed.ac.uk/handle/10283/3433, upon request to the authors.

load("/nfs/team205/ed6/data/Ramachandran2019_liver/tissue.rdata")

## Convert to SingleCellExperiment
liver_sce <- SingleCellExperiment(assay = list(counts=tissue@raw.data, logcounts=tissue@data), 
                                  colData = tissue@meta.data)

Feature selection

dec_liver <- modelGeneVar(liver_sce)

Error in modelGeneVar(liver_sce) : object 'liver_sce' not found

Dim reduction

liver_sce <- runPCA(liver_sce, subset_row=hvgs, ncomponents=11)


plotPCA(liver_sce, colour_by="condition", ncomponents=3)

scater::plotUMAP(liver_sce, colour_by="condition", point_alpha=1,  point_size=0.5) 
scater::plotUMAP(liver_sce, colour_by="dataset", point_alpha=0.3,  point_size=0.5)

scater::plotUMAP(liver_sce, colour_by="annotation_lineage", point_alpha=0.3,  point_size=0.5, text_by='annotation_lineage')

scater::plotUMAP(liver_sce, colour_by='annotation_indepth', point_alpha=0.3,  point_size=0.5, text_by='annotation_indepth')

Apply Milo

Let’s test for differential abundance between healthy and cirrhotic livers.

Sample neighbourhoods

liver_milo <- Milo(liver_sce)

## Build KNN graph
liver_milo <- buildGraph(liver_milo, d = 11, k=30)

Constructing kNN graph with k:30
Retrieving distances from 30 nearest neighbours

## Compute neighbourhoods with refined sampling
liver_milo <- makeNhoods(liver_milo, k=30, d=11, refined=TRUE)

Checking valid object

plotNhoodSizeHist(liver_milo, bins=150)

Make design matrix

liver_meta <- as.tibble(colData(liver_milo)[,c("dataset","condition")]) 

`as.tibble()` is deprecated as of tibble 2.0.0.
Please use `as_tibble()` instead.
The signature and semantics have changed, see `?as_tibble`.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

liver_meta <- distinct(liver_meta) %>%
  mutate(condition=factor(condition, levels=c("Uninjured", "Cirrhotic"))) %>%
  column_to_rownames("dataset")

Test DA

liver_milo <- countCells(liver_milo, samples = "dataset", meta.data = data.frame(colData(liver_milo)[,c("dataset","condition")]) )

Checking meta.data validity
Setting up matrix with 4973 neighbourhoods
Counting cells in neighbourhoods

milo_res <- testNhoods(liver_milo, design = ~ condition, design.df = liver_meta, fdr.weighting = "k-distance")

Performing spatial FDR correction withk-distance weighting

## Save milo object and results
saveRDS(liver_milo,"/nfs/team205/ed6/data/Ramachandran2019_liver/tissue_milo.RDS")
write_csv(milo_res, "/nfs/team205/ed6/data/Ramachandran2019_liver/milo_results.csv")

liver_milo <- readRDS("/nfs/team205/ed6/data/Ramachandran2019_liver/tissue_milo.RDS")
milo_res <- read_csv("/nfs/team205/ed6/data/Ramachandran2019_liver/milo_results.csv")

Parsed with column specification:
cols(
  logFC = col_double(),
  logCPM = col_double(),
  F = col_double(),
  PValue = col_double(),
  FDR = col_double(),
  Nhood = col_double(),
  SpatialFDR = col_double()
)

Visualize results

milo_res %>%
  ggplot(aes(logFC, -log10(SpatialFDR))) + 
  geom_point(size=0.4) +
  geom_hline(yintercept = -log10(0.1))

Check cell type composition of DA neighbourhoods

#' Add annotation of most frequent cell type per nhood to milo results table
add_nhood_coldata_to_res <- function(liver_milo, milo_res, coldata_col){
  nhood_counts <- sapply(seq_along(nhoods(liver_milo)), function(x) table(colData(liver_milo)[as.vector(nhoods(liver_milo)[[x]]), coldata_col]))
  nhood_counts <- t(nhood_counts)
  rownames(nhood_counts) <- seq_along(nhoods(liver_milo))
  max_val <- apply(nhood_counts, 1, function(x) colnames(nhood_counts)[which.max(x)])
  max_frac <- apply(nhood_counts, 1, function(x) max(x)/sum(x))
  milo_res[coldata_col] <- max_val
  milo_res[paste0(coldata_col, "_fraction")] <- max_frac
  return(milo_res)
}

milo_res <- add_nhood_coldata_to_res(liver_milo, milo_res, "annotation_lineage")
milo_res <- add_nhood_coldata_to_res(liver_milo, milo_res, "annotation_indepth")

This shows that I can decover all the clusters where DA was detected in the original paper (see all the barplots for each lineage) and more! All in a single analysis, and without knowing where the subclusters are.

milo_res %>%
  mutate(is_signif = ifelse(SpatialFDR < 0.1, 1, 0)) %>%
  mutate(logFC_color = ifelse(is_signif==1, logFC, NA)) %>%
  arrange(annotation_lineage) %>%
  mutate(Nhood=factor(Nhood, levels=unique(Nhood))) %>%
  ggplot(aes(annotation_indepth, logFC, color=logFC_color)) +
  scale_color_gradient2() +
  ggbeeswarm::geom_quasirandom(alpha=1) +
  coord_flip() +
  facet_wrap(annotation_lineage~., scales="free") +
  theme_bw(base_size=16)

Visualization with co-embedding

Instead of using the rotation matrix from original PCA, I am concatenating the nhood and sc matrix and recomputing PCA.

nhoodExpression(liver_milo)

1 x 1 sparse Matrix of class "dsCMatrix"
      
[1,] .

# pca_mat <- rbind(merged_pca$x[(ncol(liver_milo)+1):(ncol(liver_milo)+(length(nhoods(liver_milo)))),], merged_pca$x[colnames(liver_milo),])
merged_pca$x

                   PC1           PC2           PC3           PC4           PC5           PC6
    [1,] -5.746957e+00 -3.964210e+00 -2.889072e-01 -3.481046e+00 -9.266978e+00  9.852946e-01
    [2,]  2.020409e+01 -1.195912e+00  1.009486e+01  3.721923e+00 -1.159365e+00  3.910588e+00
    [3,] -8.480160e+00 -3.481662e+00  1.323102e+00  1.962629e+00  7.374387e+00  3.539573e+00
    [4,] -1.057388e+01 -3.847828e+00  3.530249e+00  7.055860e+00  1.337170e+01 -3.798475e+00
    [5,] -6.651522e+00 -3.299075e+00  8.316448e-01 -1.049599e+00 -4.691050e+00  1.607165e-01
    [6,] -2.256101e+00 -4.085494e-01 -1.735948e+00 -2.658254e+00 -5.632104e+00  3.380118e+00
    [7,] -5.286012e+00  1.419056e+00 -1.272377e+00 -3.811042e+00 -9.041807e+00  1.080641e+00
    [8,]  2.717785e+00  1.689445e+01 -4.774821e+00 -6.150700e-01  3.269642e+00  1.196702e+01
    [9,]  1.319730e+01 -9.917881e+00 -2.067499e+01  5.568129e+00 -1.504795e+00  3.685313e+00
   [10,] -8.456725e+00 -2.954216e+00  1.231968e+00 -1.438493e+00 -5.187318e+00  2.395585e-01
   [11,] -1.050696e+01 -3.671115e+00  3.057666e+00  3.357576e+00  6.949990e+00 -3.157156e+00
   [12,] -9.653671e+00 -3.806361e+00  2.157379e+00 -4.935573e-02 -2.073005e+00  9.743742e-02
   [13,]  1.341612e+01 -5.134259e+00 -7.678515e+00 -5.645511e+00 -1.785058e+00 -2.272457e+00
   [14,] -7.398999e+00 -3.492389e+00 -1.604773e-01 -3.073152e+00 -5.727850e+00  5.235551e+00
   [15,] -1.056371e+01 -3.610760e+00  2.648497e+00  2.692982e+00  5.349519e+00 -2.843876e+00
   [16,]  1.699556e+01 -3.343277e+00  1.032622e+01  1.739819e+00 -5.774999e+00 -3.081407e+00
   [17,]  2.032175e+01 -2.894313e+00  1.233561e+01  2.323507e+00 -4.797152e+00 -4.791304e+00
   [18,] -3.993270e+00  1.316032e+00 -1.632342e+00 -3.438243e+00 -9.205952e+00  7.667986e-01
   [19,] -1.017250e+01 -4.334764e+00  2.172916e+00  3.087259e-01 -1.193059e+00 -7.452379e-01
   [20,] -9.773156e+00 -1.866452e+00  2.515511e+00  1.595005e+00  1.206122e+00 -2.242694e+00
   [21,]  1.478546e+01 -7.572919e+00 -2.148589e+01  8.181546e+00 -5.286977e-01 -1.319204e+00
   [22,] -2.510969e+00  6.966758e+00 -3.273229e+00 -3.746516e+00 -5.141421e+00 -1.089334e-01
   [23,]  1.919755e+01 -3.419177e+00  1.279179e+01  2.073615e+00 -6.448676e+00 -6.272484e+00
   [24,]  7.372441e+00  2.779465e+01 -5.233624e+00  2.662509e+00  8.436619e+00  7.232789e+00
   [25,] -6.123660e-01  1.630314e+01 -4.854078e+00 -2.392557e+00 -2.487947e+00 -4.210559e+00
   [26,] -1.090812e+01 -4.592305e+00  2.418826e+00  5.432502e-01 -1.702320e+00 -3.086971e+00
   [27,]  2.384344e+01 -2.508335e+00  1.315235e+01  2.726097e+00  7.656815e-01 -7.523390e+00
   [28,]  2.128186e+01 -1.889362e+00  1.429585e+01  5.064510e+00 -3.475676e+00 -5.417601e+00
   [29,] -6.438535e+00 -3.355929e+00 -7.108547e-02 -3.265273e+00 -8.097027e+00  1.412250e+00
   [30,]  1.256496e+01 -9.657688e+00 -2.024242e+01  7.153231e+00 -2.740818e+00 -9.867003e-01
   [31,] -4.836756e+00  2.074187e+00 -1.670703e+00 -3.967458e+00 -8.979641e+00  7.138328e-01
   [32,] -5.279551e+00 -3.793559e+00 -6.009438e-01 -3.759817e+00 -9.485460e+00  1.704818e+00
   [33,] -8.390113e+00 -4.018079e+00  1.498943e+00 -1.742012e-01 -3.433592e+00 -8.341746e-01
   [34,] -8.483980e+00 -3.504860e+00  1.520554e+00 -1.247101e+00 -5.025737e+00 -8.601169e-01
   [35,] -1.524611e-01  1.329435e+01 -5.455816e+00 -3.819411e+00 -3.543739e+00  1.050834e+00
   [36,] -9.036937e+00 -3.691337e+00  1.923373e+00  3.048165e+00  5.450152e+00 -2.695680e+00
   [37,] -2.589043e-02  1.690999e+01 -4.803117e+00 -2.220435e+00 -2.333972e+00 -4.155279e+00
   [38,] -7.413219e+00 -3.309764e+00  3.135559e-01 -2.073538e+00 -2.811864e+00  6.434027e+00
   [39,]  1.947352e+01 -9.477198e-01  1.143566e+01  3.417096e+00 -3.373256e+00 -4.202561e+00
   [40,]  1.871736e+01 -8.396408e+00 -2.858435e+01  1.172955e+01  8.527806e-01  7.568445e+00
   [41,]  6.774763e+00  2.547803e+01 -4.685519e+00  2.231578e+00  7.083490e+00  4.581679e+00
   [42,]  1.921574e+01 -9.910676e+00 -2.405503e+00 -2.930826e+01  1.089380e+01 -2.713007e+00
   [43,]  2.123531e+00 -6.543081e-01 -3.309896e+00 -4.028290e+00 -3.146647e+00  1.383247e+01
   [44,]  1.974702e+01 -4.162741e+00  7.772592e+00 -2.589497e+00 -2.248232e+00 -4.813119e+00
   [45,] -1.091601e+01 -3.467889e+00  3.803579e+00  5.453100e+00  1.329413e+01 -1.344405e+00
   [46,]  3.887841e-01  1.503894e+01 -3.996655e+00 -2.030411e+00 -3.455753e+00 -5.878443e+00
   [47,] -9.546027e+00 -3.978208e+00  2.771101e+00  3.956113e+00  8.392452e+00 -1.053945e+00
   [48,]  1.904150e+01 -2.798299e+00  7.355254e+00  1.293183e+00 -3.196007e+00 -9.929789e-01
   [49,]  1.796646e+01 -2.475386e+00  1.040336e+01  1.029843e+00 -5.339073e+00 -4.431123e+00
   [50,]  1.826099e+01 -3.085680e+00  1.014878e+01  1.194658e+00 -2.592394e+00 -4.805816e+00
   [51,] -7.638750e+00 -1.591447e+00  9.420861e-01 -1.176422e+00 -4.829198e+00 -1.396014e+00
   [52,] -1.010172e+01 -2.886213e+00  2.088007e+00  6.374846e-02 -2.782631e+00 -2.172227e+00
   [53,]  2.304762e+01 -1.008835e+01 -1.477895e+00 -3.139194e+01  1.564655e+01 -3.652444e+00
   [54,] -6.014653e+00  2.259342e+00 -1.193865e+00 -2.961855e+00 -8.338538e+00 -1.564557e+00
   [55,]  2.385362e+01 -8.394732e-01  1.414020e+01  5.910616e+00  4.504022e-02  5.366587e+00
   [56,]  5.519430e+00  2.529865e+01 -5.404081e+00  1.497360e+00  6.250739e+00  5.458218e+00
   [57,] -2.272075e+00  1.533649e+01 -3.177392e+00 -1.487070e+00 -2.131017e+00 -8.274836e+00
   [58,] -8.729422e+00 -3.748911e+00  7.965584e-01 -1.303092e+00 -2.989684e+00  2.503175e+00
   [59,] -6.302169e+00 -2.971210e+00 -3.773394e-01 -3.411293e+00 -5.713570e+00  3.726058e+00
   [60,] -1.005799e+01 -4.385484e+00  1.925488e+00 -7.348042e-01 -3.793108e+00 -1.322142e+00
   [61,]  1.624498e+00  2.222640e+01 -4.611967e+00  1.018211e-01  1.812284e+00 -5.690169e+00
   [62,]  5.563327e-01 -9.008543e-01 -2.127325e+00 -3.951707e+00 -6.464204e+00  7.926228e+00
   [63,] -8.377186e+00 -2.872042e+00  1.248036e+00 -2.209730e+00 -7.892889e+00 -1.385429e+00
   [64,] -1.075776e+01 -1.918528e+00  3.444132e+00  5.007846e+00  1.063460e+01 -2.887365e+00
   [65,]  2.253735e+01 -1.511113e+00  1.504748e+01  6.092721e+00 -2.455832e+00 -3.469060e+00
   [66,] -6.927576e+00 -3.894473e+00  1.187132e-01 -1.968148e+00 -2.623587e+00  5.829123e+00
   [67,] -7.558359e+00 -3.189628e+00  3.998701e-01 -1.777296e+00 -1.929691e+00  5.414244e+00
   [68,] -9.226557e+00 -4.542261e+00  2.095193e+00  3.498706e+00  6.233660e+00 -1.607899e+00
   [69,] -7.856225e+00 -2.579283e+00  1.163552e+00 -7.917642e-01 -4.682815e+00 -7.428327e-01
   [70,] -9.308620e+00 -3.213499e+00  1.837130e+00  2.422442e+00  7.012309e+00  1.797187e+00
   [71,]  1.903127e+01 -1.259010e+01 -2.947056e+01  1.218895e+01  4.611931e-01 -8.090750e-01
   [72,] -2.913821e+00  8.373820e+00 -3.027005e+00 -2.834921e+00 -3.815701e+00 -2.343744e+00
   [73,] -8.281531e+00 -3.430443e+00  1.424895e+00  2.081903e+00  7.461858e+00  4.016385e+00
   [74,] -9.487179e+00 -4.162587e+00  1.998701e+00  5.421667e-01 -6.902656e-01 -1.830132e+00
   [75,] -9.179186e+00 -3.163205e+00  1.186173e+00 -9.854946e-01 -2.585261e+00  3.022395e+00
   [76,] -1.090814e+01 -3.026649e+00  2.515502e+00  4.896667e-01  2.534386e-01 -6.133216e-02
   [77,] -2.093092e-01  1.346187e+01 -5.639983e+00 -3.474260e+00 -3.191219e+00  3.738933e-02
   [78,]  1.293702e+01 -8.697600e+00 -1.809452e+01  4.897273e+00  9.058165e-01  6.353597e+00
   [79,] -1.078850e+01 -4.953585e+00  2.276301e+00  2.732882e-03 -2.891011e+00 -1.802549e+00
   [80,]  5.831285e+00  2.564361e+01 -5.704128e+00  1.181471e+00  7.196629e+00  8.096020e+00
   [81,] -5.938305e+00 -4.102327e+00  8.371287e-02 -2.436694e+00 -7.515218e+00  2.256787e-01
   [82,] -6.924419e+00 -3.745239e+00 -4.243890e-02 -2.023495e+00 -3.313384e+00  4.847640e+00
   [83,] -4.623047e+00  1.118727e+00 -1.804994e+00 -3.637109e+00 -9.596640e+00  1.813324e-01
   [84,] -9.211082e+00 -4.009377e+00  2.121923e+00  4.015614e+00  8.601155e+00 -1.355264e+00
   [85,] -8.766429e-01  1.837630e+01 -3.984870e+00 -1.265116e+00 -1.121295e+00 -6.776690e+00
   [86,]  2.306503e+01 -8.492689e-01  1.529742e+01  4.593190e+00 -7.322384e-01 -7.870986e+00
   [87,] -8.124183e+00 -3.453227e+00  1.132423e+00  1.322108e+00  6.019896e+00  4.460491e+00
   [88,] -5.837409e+00 -3.040563e+00 -8.976178e-01 -3.640996e+00 -7.378903e+00  4.391836e+00
   [89,] -7.165352e+00 -3.223467e+00  7.459761e-02 -2.973201e+00 -6.090416e+00  3.379939e+00
   [90,]  1.287191e+01 -1.033794e+01 -2.065140e+01  4.883429e+00 -1.542512e+00  3.027705e+00
                   PC7           PC8           PC9          PC10          PC11
    [1,]  2.859882e-01 -2.130222e+00 -5.193922e-01  2.861491e-02  4.556250e+00
    [2,] -1.137823e+00  5.726144e-01 -1.812151e+00  3.647566e-01  5.494351e-01
    [3,] -5.524589e+00  2.103488e+00  1.078860e+00  7.328789e-01 -3.430851e+00
    [4,] -1.375097e+00 -2.903202e-01 -1.458173e+00 -1.554236e+00 -2.697686e+00
    [5,]  5.749401e+00  3.002275e-01  1.182957e+00  2.510053e+00 -2.749457e+00
    [6,] -3.242642e+00 -2.552128e+00 -1.446976e+00 -8.823207e+00 -1.348790e+00
    [7,] -8.408118e-01  4.869617e-01 -1.587462e+00 -1.536191e+01 -6.630300e+00
    [8,] -2.265194e+00 -4.969788e+00  8.961229e+00  1.542750e+00  2.318473e+00
    [9,] -3.006973e+00  9.988710e-01  8.262569e-01 -4.535231e-01  9.683868e-01
   [10,]  4.714068e+00 -8.040763e-02  1.333264e+00  7.521468e-01 -1.849990e+00
   [11,] -2.216957e+00 -1.513302e+00 -2.448240e+00 -2.439649e+00  5.712344e+00
   [12,]  5.765269e+00  5.388207e-01  2.392514e+00  3.293436e+00 -2.957717e+00
   [13,]  6.608583e-01 -7.011972e+00 -5.307209e+00 -1.222657e-01  1.041111e+01
   [14,] -1.713022e+00  8.954450e-01  1.503074e+00  2.865243e+00 -2.513375e-02
   [15,] -2.075778e+00 -1.856969e+00 -2.510550e+00 -2.790289e+00  6.746634e+00
   [16,] -5.643651e+00 -3.717198e+00  1.933953e+00 -2.697930e-01  3.093029e+00
   [17,] -4.481844e+00 -1.874020e+00  2.125870e+00 -3.372444e-03  9.116394e-01
   [18,] -2.260923e+00 -1.036743e+00 -2.942329e+00 -1.381177e+01 -3.675916e+00
   [19,]  3.433053e+00 -1.495189e+00  6.209300e-01  1.197178e+00  2.112462e+00
   [20,]  4.027443e+00 -1.331553e+00  1.205616e-01 -2.697804e-01  2.231002e+00
   [21,]  4.419189e-01  1.198075e+00  4.952893e-01 -7.979612e-01  8.919654e-01
   [22,] -4.709105e+00  3.386354e+00 -1.277404e+00 -1.237353e+00 -6.432999e-01
   [23,] -4.610437e+00 -3.160656e+00  1.550091e+00 -4.635351e-01  4.255031e+00
   [24,]  6.715770e+00 -8.971997e+00  1.158366e+01 -2.552915e+00  2.649160e+00
   [25,] -3.481582e+00  5.362464e+00 -4.648477e+00  4.333401e+00 -9.727940e-01
   [26,]  5.570632e+00 -1.103419e+00  4.504034e-01  1.354254e+00  2.029751e+00
   [27,] -3.542012e+00 -1.272209e+00  8.386769e+00  4.752165e+00 -8.993496e+00
   [28,] -2.108336e+00 -2.110179e+00  2.750342e+00 -4.306310e-01 -6.813197e-01
   [29,]  1.179093e+00 -7.315015e-01  2.884892e-01 -2.452375e-02  1.481939e+00
   [30,] -6.552662e-01  2.113329e+00  2.110972e+00 -1.837788e+00  4.229650e+00
   [31,] -7.852203e-01  5.043212e-01 -1.643632e+00 -1.291813e+01 -5.357930e+00
   [32,] -4.201402e-01 -1.652320e+00 -3.477568e-01  8.171901e-02  4.225339e+00
   [33,]  5.977403e+00 -9.813005e-01  9.550793e-01  2.969373e+00 -1.115039e+00
   [34,]  6.408966e+00  2.524766e-01  1.361039e+00  2.191698e+00 -1.925151e+00
   [35,] -8.397866e+00  5.429139e+00 -6.349261e+00  8.246586e+00  5.455420e-01
   [36,] -1.818442e+00 -1.145242e+00 -1.743136e+00 -1.903365e+00  2.288128e+00
   [37,] -3.754456e+00  5.227447e+00 -5.327655e+00  5.094786e+00 -4.633104e-01
   [38,] -3.952078e+00  6.215110e-03  1.979807e+00  2.396348e+00  5.165680e-01
   [39,] -1.823137e+00 -2.074363e+00  1.739231e+00 -4.974171e-01  7.296852e-01
   [40,] -6.184990e-01 -4.840530e+01 -3.423069e+01  1.279313e+01 -2.286637e+01
   [41,]  5.887510e+00 -6.525249e+00  8.798580e+00 -1.079189e+00  1.952284e+00
   [42,]  4.413173e+00 -1.081669e+00 -8.648142e-01 -3.169551e+00  6.744029e-01
   [43,] -8.873905e+00 -3.262093e+00  3.602578e-01 -4.726203e+00 -5.501042e+00
   [44,] -3.255131e+00 -2.207620e+00 -1.286031e-01  8.718990e-01  4.618224e+00
   [45,] -1.462441e+00  1.587067e+00  7.417243e-01 -2.116672e+00 -5.657817e+00
   [46,] -2.307070e+00  2.200195e+00 -4.175007e+00  6.784984e-01  1.606286e+00
   [47,] -1.728730e+00 -3.999100e-01  3.533867e-02 -2.359918e+00 -1.179167e+00
   [48,] -1.316674e+00 -1.144927e+00 -2.137333e+00 -4.706853e-01  5.534458e+00
   [49,] -5.569837e+00 -3.081360e+00  2.879514e+00 -1.167880e-03  1.010684e+00
   [50,] -5.825921e+00 -2.562783e+00  6.552980e+00  2.788878e+00 -4.924513e+00
   [51,]  5.576843e+00 -5.483632e-01  6.125606e-01  2.762311e-01 -9.253945e-01
   [52,]  6.163084e+00 -7.615518e-01  6.233908e-01  1.087099e+00  6.288579e-02
   [53,]  6.975897e+00  1.149276e+00 -2.183843e-01 -1.172909e+00 -3.195530e+00
   [54,]  9.389802e-01  3.876052e-01 -2.467477e+00 -1.576602e+01 -5.899248e+00
   [55,]  3.290225e+00  4.339037e+00 -3.335519e+00  3.666189e-01 -2.646767e+00
   [56,]  5.431102e+00 -6.880761e+00  9.747268e+00 -2.595456e+00  2.377829e+00
   [57,]  1.233224e+00  4.359092e+00 -3.501280e+00 -4.687929e+00 -3.931292e+00
   [58,]  1.621469e+00  1.846750e+00  2.195483e+00  4.037672e+00 -1.976229e+00
   [59,] -9.508633e-02  1.567763e+00  1.703315e+00  3.941418e+00 -9.157604e-01
   [60,]  4.252696e+00 -6.747343e-01  8.011059e-01  1.159408e+00  1.759141e+00
   [61,]  2.743978e+00  2.818932e-01  1.326605e-01 -6.883191e-01 -8.753245e-01
   [62,] -6.630902e+00 -3.483951e+00 -5.227485e-01 -7.590372e+00 -1.886549e+00
   [63,]  5.624217e+00 -1.098537e+00 -9.584134e-02 -1.281254e+00  8.450182e-01
   [64,]  3.567701e-01  6.345243e-01  3.430271e-01 -1.528162e+00 -2.797286e+00
   [65,] -9.284462e-01 -8.129640e-01  8.322186e-01  4.205349e-01 -6.228501e-01
   [66,] -3.619646e+00  1.462423e-01  1.290584e+00  2.625164e+00  9.122604e-01
   [67,] -3.687939e+00  7.827507e-01  1.996782e+00  3.007815e+00 -3.839275e-01
   [68,] -1.529336e+00 -9.086508e-01 -1.011362e+00 -6.687775e-01 -1.260330e+00
   [69,]  5.997910e+00  5.430498e-02  1.056168e+00  1.943014e+00 -1.717495e+00
   [70,] -5.889380e+00  1.581000e+00 -7.204928e-02 -2.042719e-01 -5.811463e-01
   [71,]  2.030762e+00 -7.738408e+00 -5.576293e+00  1.468934e+00 -3.180082e+00
   [72,] -3.180856e+00  3.620497e+00 -1.239711e+00 -2.007913e+00 -1.348997e+00
   [73,] -6.995377e+00  1.777505e+00  1.138435e+00  4.471623e-01 -2.009539e+00
   [74,]  1.995410e+00 -1.978638e+00 -1.245585e+00 -5.425819e-03  5.657567e+00
   [75,] -1.119704e-01 -4.522286e-01  8.598876e-01  2.277517e+00  1.215323e+00
   [76,]  4.560359e+00 -4.331545e-01  2.049028e+00  2.103770e+00 -2.058371e+00
   [77,] -8.362894e+00  5.276581e+00 -8.443526e+00  9.033287e+00  1.012730e+00
   [78,] -3.963887e+00  3.840823e+00  3.211265e+00  7.140609e-01 -1.998592e+00
   [79,]  5.306629e+00 -1.226874e+00  8.618405e-01  1.267618e+00  1.969467e+00
   [80,]  2.970449e+00 -4.377120e+00  8.714987e+00  3.116594e+00  6.432366e-01
   [81,]  3.697581e+00 -1.227905e+00  2.728196e-01  5.462314e-01  9.921183e-01
   [82,] -2.213776e-01  1.604744e+00  2.660857e+00  3.877094e+00 -2.224217e+00
   [83,] -1.387621e+00 -4.704942e-01 -3.119666e+00 -1.648780e+01 -4.905789e+00
   [84,] -3.724954e+00 -2.391163e+00 -2.937785e+00 -1.426617e+00  5.991375e+00
   [85,] -1.781024e+00  4.132139e+00 -6.287664e+00  2.964170e+00  5.625553e-01
   [86,] -3.900708e+00 -2.214195e+00  8.165680e+00  2.553659e+00 -6.647968e+00
   [87,] -7.052582e+00  1.854309e+00  1.016399e+00  1.137819e+00 -1.827862e+00
   [88,] -1.374508e+00  7.317936e-01  7.304412e-01  2.312534e+00  4.470790e-01
   [89,] -2.567184e-01  1.076468e+00  9.514374e-01  3.398630e+00  3.236781e-01
   [90,] -2.182193e+00  5.674011e+00  4.166364e+00 -2.435596e-01  1.270903e+00
 [ reached getOption("max.print") -- omitted 63241 rows ]

I had to hard code the plotting function here, because using the function caused a C stack usage error, that I am yet to fix.

split_by=NULL
## Join test results and dimensionality reductions
rdim_df <- data.frame(nhoodReducedDim(liver_milo, "UMAP")[,c(1,2)])
colnames(rdim_df) <- c('X','Y')

n_nhoods <- length(nhoods(liver_milo))
rdim_df[,"Nhood"] <- ifelse(1:nrow(rdim_df) %in% c(1:n_nhoods), c(1:n_nhoods), NA)
viz_df  <- left_join(rdim_df, milo_res, by="Nhood")
viz_df[["nhIndex"]] <- unlist(ifelse(!is.na(viz_df$Nhood), nhoodIndex(liver_milo)[viz_df$Nhood],NA))
viz_df[is.na(viz_df["nhIndex"]),'nhIndex'] <- 1:ncol(liver_milo) # Add index also to single-cells

if (!is.null(split_by)){
  split_df <- data.frame(split_by=colData(liver_milo)[,split_by])
  split_df[,"nhIndex"] <- 1:nrow(split_df)
  viz_df  <- left_join(viz_df, split_df, by="nhIndex")
}

filter_alpha=0.1
## Filter significant DA nhoods
if (!is.null(filter_alpha)) {
  if (filter_alpha > 0) {
    viz_df <- mutate(viz_df, logFC = ifelse(SpatialFDR > filter_alpha, NA, logFC))
  }
}

## Plot
pt_size=1
nhood_reduced_dims="UMAP"
components=c(1,2)
  pl <-
    ggplot(data = arrange(viz_df, abs(logFC)),
           aes(X, Y)) +
    geom_point(aes(color = ''), size = pt_size / 3, alpha = 0.5) +
    geom_point(
      data = . %>% filter(!is.na(SpatialFDR)),
      aes(fill = logFC),
      size = pt_size,
      stroke = 0.1,
      # colour="black",
      shape = 21
    ) +
    scale_fill_gradient2(
      midpoint = 0,
      high = "red",
      low = "blue",
      name = "log-FC"
    ) +
    xlab(paste(nhood_reduced_dims, components[1], sep="_")) +
    ylab(paste(nhood_reduced_dims, components[2], sep="_"))

  if (!is.null(split_by)) {
    pl <- pl + facet_wrap(split_by~.)
  }
  if (!is.null(filter_alpha)) {
    pl <- pl +
      scale_color_manual(values = 'grey', label = paste("SpatialFDR >", round(filter_alpha, 2))) +
      guides(colour = guide_legend(
        '',
        override.aes = list(
          shape = 21,
          colour = "black",
          fill = "grey50",
          size = pt_size,
          alpha = 1,
          stroke = 0.1
        )
      ))
  } else {
    pl <- pl +
      scale_color_manual(values = 'grey') +
      guides(color="none")
  }

  pl <- pl +
    theme_classic(base_size = 16) +
    theme(
      axis.ticks = element_blank(),
      axis.text = element_blank(),
      plot.title = element_text(hjust = 0.5)
    )
  

pl

Running Milo one lineage at the time

Endothelial lineage

Close up in the UMAP (I don’t recompute because it changes the underlying KNN graph, and I am not interested in doing this now)

Differential expression between DA neighbourhoods

I will try this on the endothelial lineage. I want to merge overlapping nhoods with significant DA and the same direction of logFC, and then do DA between cells in up and down nhoods (I guess you could also do up or down VS all the rest). This allows us to perform a comparison without further clustering.

endo_res <- filter(milo_res, annotation_lineage=="Endothelia")

head(data.frame(markers$negLogFC_2), 20)

head(data.frame(markers$negLogFC_1), 20)

head(data.frame(markers$posLogFC_1), 20)

Visualize some of the markers that differentiate DA neighbourhoods

.plot_nhood_expression(liver_milo, endo_nhoods, features=head(rownames(markers$negLogFC_1), 40))

Joining, by = "Nhood"

Normalization

```r
lib_sf <- librarySizeFactors(liver_sce)
hist(log10(lib_sf), xlab=\Log10[Size factor]\, col='grey80')

<!-- rnb-source-end -->

<!-- rnb-plot-begin -->

<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArwAAAGwCAYAAABLkLalAAAEGWlDQ1BrQ0dDb2xvclNwYWNlR2VuZXJpY1JHQgAAOI2NVV1oHFUUPrtzZyMkzlNsNIV0qD8NJQ2TVjShtLp/3d02bpZJNtoi6GT27s6Yyc44M7v9oU9FUHwx6psUxL+3gCAo9Q/bPrQvlQol2tQgKD60+INQ6Ium65k7M5lpurHeZe58853vnnvuuWfvBei5qliWkRQBFpquLRcy4nOHj4g9K5CEh6AXBqFXUR0rXalMAjZPC3e1W99Dwntf2dXd/p+tt0YdFSBxH2Kz5qgLiI8B8KdVy3YBevqRHz/qWh72Yui3MUDEL3q44WPXw3M+fo1pZuQs4tOIBVVTaoiXEI/MxfhGDPsxsNZfoE1q66ro5aJim3XdoLFw72H+n23BaIXzbcOnz5mfPoTvYVz7KzUl5+FRxEuqkp9G/Ajia219thzg25abkRE/BpDc3pqvphHvRFys2weqvp+krbWKIX7nhDbzLOItiM8358pTwdirqpPFnMF2xLc1WvLyOwTAibpbmvHHcvttU57y5+XqNZrLe3lE/Pq8eUj2fXKfOe3pfOjzhJYtB/yll5SDFcSDiH+hRkH25+L+sdxKEAMZahrlSX8ukqMOWy/jXW2m6M9LDBc31B9LFuv6gVKg/0Szi3KAr1kGq1GMjU/aLbnq6/lRxc4XfJ98hTargX++DbMJBSiYMIe9Ck1YAxFkKEAG3xbYaKmDDgYyFK0UGYpfoWYXG+fAPPI6tJnNwb7ClP7IyF+D+bjOtCpkhz6CFrIa/I6sFtNl8auFXGMTP34sNwI/JhkgEtmDz14ySfaRcTIBInmKPE32kxyyE2Tv+thKbEVePDfW/byMM1Kmm0XdObS7oGD/MypMXFPXrCwOtoYjyyn7BV29/MZfsVzpLDdRtuIZnbpXzvlf+ev8MvYr/Gqk4H/kV/G3csdazLuyTMPsbFhzd1UabQbjFvDRmcWJxR3zcfHkVw9GfpbJmeev9F08WW8uDkaslwX6avlWGU6NRKz0g/SHtCy9J30o/ca9zX3Kfc19zn3BXQKRO8ud477hLnAfc1/G9mrzGlrfexZ5GLdn6ZZrrEohI2wVHhZywjbhUWEy8icMCGNCUdiBlq3r+xafL549HQ5jH+an+1y+LlYBifuxAvRN/lVVVOlwlCkdVm9NOL5BE4wkQ2SMlDZU97hX86EilU/lUmkQUztTE6mx1EEPh7OmdqBtAvv8HdWpbrJS6tJj3n0CWdM6busNzRV3S9KTYhqvNiqWmuroiKgYhshMjmhTh9ptWhsF7970j/SbMrsPE1suR5z7DMC+P/Hs+y7ijrQAlhyAgccjbhjPygfeBTjzhNqy28EdkUh8C+DU9+z2v/oyeH791OncxHOs5y2AtTc7nb/f73TWPkD/qwBnjX8BoJ98VQNcC+8AAEAASURBVHgB7d0HmFPF2sDxd4FlAem9KbCgiKiIigqiiKKCBVBZ9QICFsSrYgPFhggK3isqCl4/e0FRUbAiogKCXqxIE1BBei/S+wLnO+9834nJbnaT7OZszkn+8zxLklPmzPwmu7yZzMxJs+wkJAQQQAABBBBAAAEEklSgWJLWi2ohgAACCCCAAAIIIGAECHh5IyCAAAIIIIAAAggktQABb1I3L5VDAAEEEEAAAQQQIODlPYAAAggggAACCCCQ1AIEvEndvFQOAQQQQAABBBBAgICX9wACCCCAAAIIIIBAUgsQ8CZ181I5BBBAAAEEEEAAAQJe3gMIIIAAAggggAACSS1AwJvUzUvlEEAAAQQQQAABBAh4eQ8ggAACCCCAAAIIJLUAAW9SNy+VQwABBBBAAAEEECDg5T2AAAIIIIAAAgggkNQCBLxJ3bxUDoG/BbZu3Srnnnuu+Rk0aNDfO4Ke/eMf/zD7r7zyysDWDRs2BM577bXXAttjfXL48GGZN29erKdxfJCAGt5///3SrFkzqVKlipx++ukyefLkoCP+fhqvdvs7x9ieZWdny2WXXSannnqqjBgxIuzJy5Ytk/79+0u7du1E33svv/yyWJYV9ljdeNZZZ8kRRxwhWjdNjz/+eOC9qTZ5bZs/f37guAkTJpjjEvFPuPZ75ZVXpGTJkqb+iSgT10QgZQTsPy4kBBBIAYH169drJGF+7EAkbI2POuoos79OnTqB/cuXLw+c9+CDDwa2x/JkxowZ1imnnGJ169YtltM4NofAo48+GmgLpy2//fbbHEf938t4tFvYjKPcePvttwfKeuedd+Y664cffrAqV64cOMapzxVXXGHZgWGu4z/++GNzbFZWVmBfr169AucfPHjQbA+37b///W/guBdeeCFwflE/yav9WrdubaWlpVlz5swp6iJxPQRSRqBEykT2VBQBBAokUKZMGenYsaM5t0mTJgXK48wzzzTnHXvssQU6n5P+T8AOEs2TEiVKyMyZM6VixYpSu3ZtT/EcOHBAnnzySXnmmWfyLde9994rW7ZskXLlykm/fv1E6zZp0iQZP368fPHFF9K+ffvA+doz+sADD5jXN910U2B7uCfNmzc3+eo+O4gMd0jCtuXVflonOyiXhx56SOzAPmHl48IIJLMAAW8yty51QyAOAtWqVcv3P+Ht27eL3XsslSpVkqpVq0qxYgUfKbVp0yZzvn5dn1+ye/Nk9erVUq9evZiCGj2nbt26IVnv3LlT1q1bJ+np6VKzZk0pXbp0yP5wL/bu3Ss6RCRnsLljxw7RgE8dCpIi1X/NmjUm22OOOcYMayjINYLPOXTokKxcuVLsnn0pXrx48K5cz6Mx14DummuukT///DPX+cEbZs+eLdOmTTObdEiDBnq7du2S+vXry19//SVPPfVUSMD7zjvviA5L0PfiOeecE5xVrue33Xab6E+ktH//fvO+1brHIzBWH20fNdX3kX5QzJnyaj/9QKnDGj755BP5+eefpUWLFjlP5TUCCBRSoOD/MxXywpyOAAL+EFi1apWULVvW/DzyyCOBQr/77rvSsGFD08uoPbc1atQwAeDIkSMDx/zyyy/mPGfD2LFjzetnn33W2SQaPGqAor3H1atXN8GiBqUa9DjjMp2D9dh//vOfplewQYMGZhzrY489Zo51yugc+/333wfK/dlnn0njxo3lyCOPNMGd9iz+/vvvYn+VLOXLlzf7MjMzTV3sr8Rl3759TjYmCHHy1qCrc+fOUqFCBbGHfYieowHK2rVrzRhR+yt6U4emTZvKr7/+GsgjvyfR1H/KlCmmLnPnzjVZadm1TG3atMkv6zz3aU9qhw4dzIcUrYOOiz3//PNl0aJFuc6JxVzLqcGuBq4vvfRSrrycDd99953z1Izf1Rdan7PPPttst4fABPbrk1dffdW8vuSSSyJ+oNLeUqe9NPjMmXbv3i09evQw7a7l1DbT95r9vW7OQ6N6reOJNT/9oKT56e+EXl899XdHU6T2017utm3bmmNHjRplHvkHAQTiLJAygzeoKAIpLhA8hve0006z3nvvvVw/ds+kGesYaQyvPVHKjDm0/xxZdq+bZQ9ZsOyv1825OhZx6tSpRvunn34y2/S44B87wDD77WDKsicshewLPs4OcEJaTV87++2eZMsOLMxrLa+z3TkheNymHYwH9tuTqCw7ELLUwDnH7qm1ND/ntR3AONlYH3zwQWB7rVq1zHHB+dk925ZeX+ttB8KBY+2g3dqzZ08gn3BPoq2//RV/IF+njPpo9wSGy9Zsy2sM78SJEy27NzFsfnbQZn3zzTchecZiru06ePBgS+u1efPmwDVyjuG1J94F9q1YsSJwvVtuuSWwfdu2bWa73fMbKO/TTz8dOFafhBuvG25b8HshIyPDXMP+oBO4lloOGTIkJO9oX9hBusnHHmZinXDCCeZH3wuapx3EmvHI0bTfPffcY87RPEgIIBB/Af1US0IAgRQQCA54g4OmcM8jBbw33nij+c9ZA0VngpHdK2rZXzdbds+hZffgGlHdtnjxYnOsXsf+6ta8tocDmP0333xzYJ8Gmb/99pulk7Ds1SQC2+2eQnOsBtFOWTXAtns5LXsIgaUT6Zzt+uik4CBHt2uwpBOfPv/8c2v69OmW3aNs2cMYAsG5PTTDcibt2T3BTjYhAa/W1x4CYPbZPc2B6+oHhSVLlphAz15tILDd7uUN5BPuSbT11wBSr3vccceZvO0edfNa2zSvFC7g1TxKlSpl8tA62qtuWHYPpfXvf/878OHh6KOPDgTqsZoHlyW/gDc4KNXrO8ke3hCw0/eCJrt3NLBNP3wEp+B8op20pm1ufzth6fFffvmlpQ76/rB7uY1FcP6Rnts9+4GyBU+Ge/vtty17LLEJyJcuXWreF5Haz+7ZNXnphxF7hYtIl2Y/AgjEKPD3/w4xnsjhCCDgL4GcAa/2aOb8cQLHSAHv3XffHfiPXldf0F49DTCdoCOnjJNvzlUanADTHrMbcq4Gj845Z5xxhskuOLC1x3QGLmGPmQ0Ea3qOk4ID3osvvtjZHPKo52qyx3Oa8p988snmuvZXzIHjgnt4hw8fHthuD88IlFF7LJ00ZsyYwPZPP/3U2Rz2MZb6awYnnXSSyTuaXsBwAa89NCBQNnupsJAy6QcOx9we8mD2xWoenGF+Aa+97F3gWvaY3cBpTi+nlsMeDmO2v/nmm4Fj7eEjgWP1SUECXv1QFpz0g5BTb3vJsuBdEZ/rBy577LM5X98zavjGG29Y9pjwsOfm134fffRRoBwLFy4Mez4bEUCg4AKM4bX/0pEQSDUBXR9Vxzfm/NEJPNEku4fXTCDSY3Wcrq7rq+NhdQxu3759RSdfRUq6Bqvd62UOswPSkElTOq7UDurMPh03a/+JEzsIDmSp67Y6SSebOeMfnW05H8NNAtIJRDp2U8fB6jheLf+sWbPMqXq9cCl4wlvw5DYdG+wkHQ/rJJ0YlVeKtf555RPLdmeimJ7TqVOnkFODXzvjjwtjHpJ5jhf6PnFSsFHw2Gmd+KXJWXNXnwc76+uCpAsuuCDkNLvXPPDa7o0NPI/miY691fe7Jp38OHr0aOnZs6cZy65rJOsktGhTcN30PU9CAIH4ChDwxteT3BBICYFGjRqZZbHsnl4zUceptE4G0wlpOmFHZ63nl4IDRl1mK2dyZs7nfNTjch6vM/vzS7p8V3DSgMLumRZdGktXFjjvvPPMzRG03Jrsnu/gwwPP7eEAgedOuXSDTnyKNcVa/1jzD3d88DVzrsoQXB/n3OBtsZo7eYR7dIJZ3WeP1Q0c4jzXsukkSE3BZbbH3waOLeiT4DbUPILbWj88xZr0Q9P7779v3vNO+fQDkz1+3XyoGDduXFRZBq/qoCs2kBBAIL4C4f+qx/ca5IYAAkkooD2ZOiNeeyq1Z+z55583qx1oVXU1AXsMZq5aB/ecatCjs9o12V/9h6zIYE9kCtyVTVc80MBLe32dZA9XcJ6K/dW5CVoDG8I8yRlA6Bqx2nOoAYr99bHoKg533HFHYCmp4CAoTHa5NgUHhrl25rEh1vrnkU1Mm1u2bBk4Pud6r8G9kU7vemHMAxcK88Rpd93l9Cbr8wULFuhDyDJpwT2wwb295sAC/KPfSAQn55q6Lbi+wcfk91y/JVEve0iCCd51RQZdScRJuppJNEmX93OSLjtHQgCB+AoQ8MbXk9wQSAkBe6yrWT5MH3UNW10irE+fPqJDE5yk6/I6yek50+W79Gtr7QnW5CxDpUMg7Bn6ooGurtFqT+ZyTpXrr7/ePNdbzzqBZe/evUWXb9IgQ5eqyrl8WeDk/3+SszfTWX5Ltzvl1K/vv/76a3OGLl1VFCmW+sejPDpswwnm7bGrpmdSv4r/z3/+I/ZYZXMJ7b13ylUY8/zKq0u7Ob3u9rho83748MMPzfAYPe+GG24InK4feJwUj4DXXp0kUFcNTp944gmTvZbHXrnDuVRUj3qTDO2B1mX57HHsor3HevtuXb7P+ZDl1DNShvq7oUnfk7q0GQkBBOIsYPe4kBBAIAUEgietFfbWwjqRyA4+zSQbnVV+1llnWXbQG5h0Ywc0IaK6qoD9p8v86HnOLYp1tQZ7rGNgn/2ffeC5Hn/RRReF5HPXXXeF7NdjdHa9rrjg5O+cEDxpLXgGve63b3QQON4eT2rZQxrMig1OnTQvZyWJ4ElrdlDmZG/phDTnmrrEm5OCJx/ZX2c7m8M+xlr//CY95bxAuElreoy2XbBz8HNdsktXsAhOsZgHn5ffpDU9buDAgQG/YHddLix4Ipse6yyX98orr+jLQCrIpDX7ZiXmusFLyGk72oFvIN9on+jqGTqB0HkfaN76u2AHwWab/m78+OOPgezya7+HH37YnGP3MgeO5wkCCMRPgB5e+y8VCQEEYhPo3r276E0k9GYRemcxeykxM7RBe7jsNVfN5J3gHO01TkPGYtrr05rd2vv11VdfybXXXmvugKZfD2vSXtdhw4blmvSjt6zVnkmdhKbn6phbe+kssdfWNec5vWrmRT7/2KsqSNeuXc0RGzduFJ3MZa8cYHqMndPs9Wqdp649xlr/eBRE207HlbZq1cr0Qqq59io6k/ac3l3nWvEyd/JzHvU9MXToUHOTBvu/NLNZe3P1phM5x0TreGtN2k6FTTqxTO8GZwerJisdmqPjcPX2xrEmfb/r+0/fvzo8Rr+h0N8FzdsObkXfQ9H2GjtDLXL6x1omjkcAgfACaRo7h9/FVgQQQCCygAaM+nWsBm+6ikHOyU1ODhoE6F247F4wsyqCsz340e6VNF+5h1stQld00IBCz7eXTQtZ1eHSSy8Ve0kps0pELF9769AKvd2wvfZsSEAeXKaifJ5f/d0ohw4v0TbRr9CDJ4c513LD3MnbedSA+48//jB3zXMmqjn7nEf7xg3mVsP6/nLuXubsK+ijfujS95OOl9WAv7BJLfW9pPnq+1d/H6JN+qFRb6ett1fWYTX20mnRnspxCCAQpQABb5RQHIYAAokV0J40XU1Bk44z1SBIl7fSCWxZWVkmWNDxvDoBjhQfAS+Za6/+zJkzzWRGZ1JdfGqZ+FycgF7f1zq+3B7ikfhCUQIEkkyAgDfJGpTqIJCsAtqDpmubzps3L2wVtZfOvjGB2He4CrufjbELeMlcV5XQyW633Xab6CobbiRdWSTSms7OdTUAj9eHK10DWVfJ0OXNunTp4lyCRwQQiKMAAW8cMckKAQTcFdDhEzqrX2f06zAKvWmBDm9o1qyZufmFM5bX3VKkVu5eMdfRd9rOOuxDhzXYk87i3hA6vENXpogmnXjiiWJPoovm0HyP0WX9tGdXl4wLXm4v35PYiQACMQsQ8MZMxgkIIOAVAWfClVfKkwrlSKR5dna26I9OFnOWV/O7uXrqBzedcJnX+He/15HyI+AFAQJeL7QCZUAAAQQQQAABBBBwTYBlyVyjJWMEEEAAAQQQQAABLwgQ8HqhFSgDAggggAACCCCAgGsCBLyu0ZIxAggggAACCCCAgBcECHi90AqUAQEEEEAAAQQQQMA1AQJe12jJGAEEEEAAAQQQQMALAgS8XmgFyoAAAggggAACCCDgmgABr2u0ZIwAAggggAACCCDgBQECXi+0AmVAAAEEEEAAAQQQcE2AgNc1WjJGAAEEEEAAAQQQ8IIAAa8XWoEyIIAAAggggAACCLgmQMDrGi0ZI4AAAggggAACCHhBgIDXC61AGRBAAAEEEEAAAQRcEyDgdY2WjBFAAAEEEEAAAQS8IEDA64VWoAwIIIAAAggggAACrgkQ8LpGS8YIIIAAAggggAACXhAg4PVCK1AGBBBAAAEEEEAAAdcECHhdoyVjBBBAAAEEEEAAAS8IlPBCISgDAgh4T2DKlCnSrl07qVKlStwLt2PHDuncubO89957cc+bDBFAAAEEEMgpkGbZKedGXiOAAAL9+vWT7du3yyWXXBJ3jK1bt8ojjzwiS5cujXveZIgAAggggEBOAXp4c4rwGgEEAgLp6elSuXLlwOt4Pdm5c6cUK8aIqnh5kg8CCCCAQP4C/I+Tvw97EUAAAQQQQAABBHwuQMDr8wak+AgggAACCCCAAAL5CxDw5u/DXgQQQAABBBBAAAGfCxDw+rwBKT4CCCCAAAIIIIBA/gIEvPn7sBcBBBBAAAEEEEDA5wIEvD5vQIqPAAIIIIAAAgggkL8AAW/+PuxFAAEEEEAAAQQQ8LkAAa/PG5DiI4AAAggggAACCOQvQMCbvw97EUAAAQQQQAABBHwu4PuAV++MvGnTJtmyZYvPm4LiI4AAAggggAACCLgh4MuAd82aNTJgwACpX7++lCxZUqpXry5VqlSRChUqSLNmzaRfv36ya9cuN7zIEwEEEEAAAQQQQMBnAiV8Vl5ZsWKFtG7dWtLS0iQrK0syMzOlcuXK5rX28i5btkzGjRsn48ePlylTpkjDhg39VkXKi0BKCBw8eFAWLFjgSl3Lli0r9erVcyVvMkUAAQQQ8J+A7wLe4cOHm57dyZMnS0ZGRljxYcOGSYcOHWT06NEyePDgsMewEQEEEiewdetW8+G1U6dOrhRiyZIl8vvvv0vjxo1dyZ9MEUAAAQT8JeC7gHfOnDnSs2fPPINd5U9PT5devXrJqFGjCHj99X6ktDEKHDhwQPbu3RvjWdEdrnm7lbTM1apVk1dffdWVS9x2222yefNmAl5XdMkUAQQQ8J+A7wLeVq1ayYwZM6R37975ak+dOlXq1KmT7zHsRMDvAjVr1pSdO3dKqVKl4l4VDUoj/Z7F/aJkiAACCCCAgAsCvgt4u3btKhr0btiwQbp162bG6OqEtWLFipmVGpYvXy5jxoyRiRMnig57ICGQzAKHDh2Sjz76SMqVKxf3at50001xz5MMEUAAAQQQSISA7wLek046SebNmyd9+vQxQxsOHz6cy61du3by5ZdfSps2bXLtYwMCCCCAAAIIIIBAagn4LuDV5mnUqJFZgUHHGOqqDXPnzhUNfHWCSt26dc0SZanVjNQWAQQQQAABBBBAIC8BXwa8ug7vyJEjZezYsaLPdXkjTeXLlzcrOGgPr67OoEsTkRBAAAEEEEAAAQRSW8B3AS/r8Kb2G5baI4AAAggggAACsQr4LuBlHd5Ym5jjEUAAAQQQQACB1Bbw3a2FdR3eHj16RLUO76RJk1K7dak9AggggAACCCCAgPgu4HXW4Y3UdqzDG0mI/QgggAACCCCAQGoI+G5IA+vwpsYbk1oigAACCCCAAALxEvBdwMs6vPFqevJBAAEEEEAAAQRSQ8B3Aa82S/A6vCtXrhS9u1p2drbUrl270Ovw/vzzz2a5s0jNP2vWLBk6dKi0bNky0qHsRwABBBBAAAEEEEiggC8DXserZMmS5tbCFSpUkOLFi0vlypWdXQV+LFOmjNSqVSvi+UuWLJG//vor4nEcgAACCCCAAAIIIJBYAV8GvG7eeKJp06aiP5HSuHHjpGrVqpEOYz8CCCCAAAIIIIBAggV8F/By44kEv2O4PAIIIIAAAggg4DMB3wW83HjCZ+8wiosAAggggAACCCRYwHfr8HLjiQS/Y7g8AggggAACCCDgMwHfBbzceMJn7zCKiwACCCCAAAIIJFjAd0MauPFEgt8xXB4BBBBAAAEEEPCZgO8CXm484bN3GMVFAAEEEEAAAQQSLOC7gFe93LzxRILbg8sjgAACCCCAAAIIxFnAlwGvY6A3ntDgV380HT58WJYtWyZ6I4oSJXxdNaeKPCKAAAIIIIAAAggUUsB3k9a0vvPnz5ebbrpJrrvuOpk+fbohePLJJ6VmzZom+K1UqZK8+OKLhaThdAQQQAABBBBAAIFkEPBdN6gGuy1atBDt3dU7nb377rvyzDPPyMMPPyxXX321nHfeefL+++/LP//5T2nQoIGcf/75ydBO1AEBBBBAAAEEEECggAK+6+H917/+Jc2bN5e1a9fKn3/+KbfccovceOON0q9fP3nppZdM0Dt+/Hhp3769PPfccwVk4TQEEEAAAQQQQACBZBHwXcC7aNEi0aXJjjjiCElLS5NrrrnGtEVWVlZIm3Tp0kWWLl0aso0XCCCAAAIIIIAAAqkn4LuAt3bt2jJ16tRASznPp02bFtimT3ToQ926dUO28QIBBBBAAAEEEEAg9QR8N4ZXJ6t16NDBjOPVMbxTpkyRu+66S4YOHSoHDx6Udu3ayeeff27G9b722mup16LUGAEEEEAAAQQQQCBEwHcBr47Nfe+990xAu337drMaQ69evWTDhg1y5513imVZZqiDPneGO4TUmBcIIIAAAggggAACKSXgu4BXW0fH6+Ycs/vWW2/J8OHDZdasWXL88cdLvXr1UqohqSwCCCCAAAIIIIBAeAFfBrzhqyJSq1Ytufjii/PazXYEEEAAAQQQQACBFBTw3aS1FGwjqowAAggggAACCCBQCAEC3kLgcSoCCCCAAAIIIICA9wUIeL3fRpQQAQQQQAABBBBAoBACBLyFwONUBBBAAAEEEEAAAe8LEPB6v40oIQIIIIAAAggggEAhBAh4C4HHqQgggAACCCCAAALeFyDg9X4bUUIEEEAAAQQQQACBQggQ8BYCj1MRQAABBBBAAAEEvC9AwOv9NqKECCCAAAIIIIAAAoUQIOAtBB6nIoAAAggggAACCHhfgIDX+21ECRFAAAEEEEAAAQQKIUDAWwg8TkUAAQQQQAABBBDwvgABr/fbiBIigAACCCCAAAIIFEKAgLcQeJyKAAIIIIAAAggg4H0BAl7vtxElRAABBBBAAAEEECiEAAFvIfA4FQEEEEAAAQQQQMD7AgS83m8jSogAAggggAACCCBQCAEC3kLgcSoCCCCAAAIIIICA9wUIeL3fRpQQAQQQQAABBBBAoBACBLyFwONUBBBAAAEEEEAAAe8LEPB6v40oIQIIIIAAAggggEAhBAh4C4HHqQgggAACCCCAAALeFyDg9X4bUUIEEEAAAQQQQACBQggQ8BYCj1MRQAABBBBAAAEEvC9AwOv9NqKECCCAAAIIIIAAAoUQ8H3Aa1mWbNq0SbZs2VIIBk5FAAEEEEAAAQQQSFYBXwa8a9askQEDBkj9+vWlZMmSUr16dalSpYpUqFBBmjVrJv369ZNdu3Yla5tRLwQQQAABBBBAAIEYBErEcKwnDl2xYoW0bt1a0tLSJCsrSzIzM6Vy5crmtfbyLlu2TMaNGyfjx4+XKVOmSMOGDT1RbgqBAAIIIIAAAgggkBgB3wW8w4cPNz27kydPloyMjLBqw4YNkw4dOsjo0aNl8ODBYY9hIwIIIIAAAggggEBqCPhuSMOcOXOkR48eeQa72mzp6enSq1cvmTRpUmq0IrVEAAEEEEAAAQQQyFPAdwFvq1atZMaMGXlWyNkxdepUqVOnjvOSRwQQQAABBBBAAIEUFfDdkIauXbuKBr0bNmyQbt26mTG6OmGtWLFiZqWG5cuXy5gxY2TixImiwx5ICCCAAAIIIIAAAqkt4LuA96STTpJ58+ZJnz59pGfPnnL48OFcLdiuXTv58ssvpU2bNrn2sQEBBBBAAAEEEEAgtQR8F/Bq8zRq1MiswHDgwAFZuXKlaK9udna21K5dW+rWrWuWKEutZqS2CCCQU0DX6NYfN5KuEkNCAAEEEPCPgC8DXodX1+DVZcd0/d3ixYub5cmcfTwigEDqCqxfv17OOuss1wCefvppuf32213Ln4wRQAABBOIr4MuAV288MXLkSBk7dqzo84MHDxqV8uXLmyXLdEiDLkdWtmzZ+GqRGwII+EJA/yY8++yzcsIJJ8S9vDo/YObMmXHPlwwRQAABBNwT8F3Ay40n3HszkDMCCCCAAAIIIJCMAr4LeLnxRDK+DakTAggggAACCCDgnoDv1uHlxhPuvRnIGQEEEEAAAQQQSEYB3wW83HgiGd+G1AkBBBBAAAEEEHBPwHdDGrjxhHtvBnJGAAEEEEAAAQSSUcB3Aa/bN56YO3eufPzxxxHbevXq1bJx48aIx3EAAggggAACCCCAQGIFfBfwKlfwjSd01QYNUvWOa40bNy70jSc0n0OHDkVsFbcWtI94YQ5AAAEEEEAAAQQQiEnAlwGvm+vwNm/eXPQnUtJbF1evXj3SYexHAAEEEEAAAQQQSLCA7wJe1uFN8DuGyyOAAAIIIIAAAj4T8F3Ayzq8PnuHUVwEEEAAAQQQQCDBAr5blox1eBP8juHyCCCAAAIIIICAzwR8F/CyDq/P3mEUFwEEEEAAAQQQSLCA74Y0sA5vgt8xXB4BBBBAAAEEEPCZgO8CXrfX4fVZ+1FcBBBAAAEEEEAAgQgCvgt4tT7B6/CuXLlSli9fLtnZ2VK7du1Cr8MbwYvdCCCAAAIIIIAAAj4T8GXA6xiXLFnSBL8aAGvSO59lZGQ4u3lEAAEEEEAAAQQQQEAiTlr79ttv5c477zR3M/OK108//SSXXHKJ7Ny50xTp008/lfr160uNGjWkYsWKcsopp4iWm4QAAggggAACCCCAQMSAVwPIL774QnTsrP48/fTTsmnTpoTJ/fjjj3LGGWcEbv/7ww8/yGWXXSZVqlQRXaN3xIgRUrZsWbngggsIehPWSlwYAQQQQAABBBDwjkDEgPeEE06QhQsXys8//yxnn322DB06VOrUqSOdO3eWDz/80IydLcrqvPPOO6Z39/PPP5dy5crJiy++KDVr1hQNhPv37y99+/aV6dOny5lnnimvv/56URaNayGAAAIIIIAAAgh4UCBiwOuU+dRTT5WRI0fK2rVr5YMPPjA9qr169TITxe644w75/fffnUNdfZw5c6a0b98+cI1t27ZJp06dpESJ0OHIunyZ3qSChAACCCCAAAIIIJDaAlEHvA6Troowe/ZsmTVrluzYsUMyMzNN7+pxxx0nAwcOdA5z7fHEE0+U999/X/bt22eucc4558gnn3wSGOKgGy3LEu0BbtKkiWvlIGMEEEAAAQQQQAABfwiEdovmUebNmzfL2LFj5a233hIdM6uTw6655hrzumnTpuYsDUKvvPJKM5725JNPziOnwm8eMGCAmZR2+umny3333Sd657Vjjz1W2rRpI9dee60Z5jBmzBgz7njGjBmFvyA5IIAAAggggAACCPhaIGLA+80330i7du1MJXVlBO1N7dChQ64hBBdddJE5ZsOGDa6C1KtXT7RMjz76qHTv3j2kZ9cJcJs3b27Kqas1kBBAAAEEEEAAAQRSWyBiwFumTBl5/PHHpVu3blKtWrU8tXT92/Xr15ve3zwPitMOHT7x9ttvy6hRo2TZsmWyZs0a2b59uxlPfOSRR0rjxo3jdCWyQQABBBBAAAEEEPC7QMSAVyeraU+pro6gy4HpmF1NN998s9x2221mOIG+1kljOtShKJMuRaY/WkYSAggggAACCCCAAALhBKKatHbhhReKrsigPbiadFLY3Llz5fjjjzfr8obLmG0IIIAAAggggAACCHhBIGLAu3jxYrOu7fz5880EMS10Wlqa6HjZp556Sh544AE5cOCAF+pCGRBAAAEEEEAAAQQQyCUQMeCdMmWKGTJwzDHH5DpZJ43t2bNHdKkyEgIIIIAAAggggAACXhSIOIZXl/z66aefZNWqVaITwoLThAkTpHjx4ubOa8HbeY4AAn8L6O/O/v37/94Qx2c6vIiEAAIIIIAAAvkLRAx4daKa3rpX19196KGHRHt69aYPeqvhBx98ULKysqR06dL5X4W9CKSogN4ZsEWLFrk+LMaLw7kBS7zyIx8EEEAAAQSSUSBiwFuqVCn57rvv5OKLL5bzzjsvxKBLly7y3HPPhWzjBQII/C2gw330boCDBw/+e2McnzlrZMcxS7JCAAEEEEAg6QQiBrxaYx3KMG/ePFm9erXMmTPHDGPQoQ4NGjRIOhAqhAACCCCAAAIIIJBcAlEFvNrD+/TTT5uAN9yKDPq1LQkBBBBAAAEEEEAAAS8KRAx4dcJN+/btpUKFCtKyZUspX768F+tBmRBAAAEEEEAAAQQQCCsQMeDVZcmKFStmhjRUqlQpbCZsRAABBBBAAAEEEEDAqwIR1+HVZY/q1q0rBLtebULKhQACCCCAAAIIIJCfQMSAt23btrJ06VJZsGBBfvmwDwEEEEAAAQQQQAABTwpEHNKgy5J17dpV2rRpI1dddZXp7dWbTQSne+65J/glzxFAAAEEEEAAAQQQ8IxAxIB37ty58v7775sCv/XWW2ELTsAbloWNCCCAAAIIIIAAAh4QiBjwXnjhhbJ9+3YPFJUiIIAAAggggAACCCAQu0DEMbzBWe7evdus1rBr1y7Zv39/8C6eI4AAAggggAACCCDgSYGoAl69w5reRrhs2bLSrFkzWbhwoQwYMED69+8ve/bs8WTFKBQCCCCAAAIIIIAAAioQMeDVO6t17NhRFi1aJCNGjJAyZcoYubPOOktefvllueOOO5BEAAEEEEAAAQQQQMCzAhHH8E6ePFnWrFljAl6929qgQYNMZa644gpz17WePXuKrtWblpbm2UpSMAQQQAABBBBAAIHUFYjYw7t48WJp2rSpubVwTqYWLVrIunXrZPny5Tl38RoBBBBAAAEEEEAAAU8IRAx4GzRoIDNmzJBNmzblKvC7774rJUqUkDp16uTaxwYEEEAAAQQQQAABBLwgEHFIw7nnnmtuNnHRRRfJXXfdJYcPHxbt9f3000/l+eefNzelKFmypBfqQhkQQAABBBBAAAEEEMglEDHg1ZUZPvzwQ+nVq5cJbjWH7t27m4w6d+4szzzzTK5M2YAAAggggAACCCCAgFcEIga8WtATTzxRZs6cKb/88ovp3dUeXR3X26RJE6/Ug3IggAACCCCAAAIIIBBWIKqAV88sVqyY6CQ1/SEhgAACCCCAAAIIIOAXgYgB7w8//CD33ntvvvWZNm1avvvd3KlLom3evFmKFy8ulStXdvNS5I0AAggggAACCCDgQ4GIAW9GRobUqlUrpGo7duyQ33//XfQObIm48YSuCzxy5EgZO3asWSP44MGDpnzly5eX+vXrS7t27WTw4MHmznAhBecFAgggEAcBvcPkihUr4pBT7iwqVqwYdhnI3EeyBQEEEEAgWoGIAW/z5s3lnXfeyZWf9qz27dtXVq5cmWufmxv0P5nWrVubG11kZWVJZmam6dnVG19s2bJFli1bJuPGjZPx48fLlClTpGHDhm4Wh7wRQCDFBPSD/gcffCDff/+9KzXXtc11NRxu5uMKL5kigECKCkQMePNy0T/G/fv3l0aNGslLL71UZL2pw4cPN724egc47X0Ol4YNGyYdOnSQ0aNHm57ecMewDQEEECiIgH6w1m+RBg4cWJDTI57Tpk0b7l4ZUYkDEEAAgdgEIt54Ir/sNmzYIIcOHZLt27fnd1hc982ZM0d69OiRZ7CrF0tPTzfLqE2aNCmu1yYzBBBAAAEEEEAAAf8JROzh1dsGf/LJJyE10yBXezm0Z/fYY48t0juttWrVytz5rXfv3iFlyvli6tSpRVqunNfnNQIIIIAAAggggIA3BCIGvAsXLpS77747V2mPOOIIOemkk4r8xhNdu3YVDXq1d7lbt25mjG6VKlXMsmkahGuAPmbMGJk4caLosAcSAggggAACCCCAQGoLRAx49ZbC+/fv94ySBtnz5s2TPn36SM+ePc3kjpyF0/F1X375pehYOBICCCCAAAIIIIBAagtEDHi9yKMT5XQFhgMHDphVIrRXNzs7W2rXri1169YV7fElIYAAAggggAACCCCgAhED3mhuPBFM+eqrr5qlwoK3ufVcb3Gsy45VqFCBG0+4hUy+CCCAAAIIIICAzwUirtJQqVIlKV26tHzzzTdm9YNmzZpJtWrVZO7cuTJ9+nQzdlZvTOH8lCgRMYYuNJneeGLAgAFmeTINeqtXr256dTXw1fL169dPdu3aVejrkAECCCCAAAIIIICA/wUiRqdly5YVXQpMJ4Cde+65gRrrSg233nqr/Prrr2FvTBE4MM5PuPFEnEHJDgEEEEAAAQQQSHKBiAGvTv7Su60FB7tqUrx4cXnwwQfNmFm9M1DO2w+75caNJ9ySJV8EEEAAAQQQQCA5BSIOaShfvrz89NNPYYcIrF+/3gxp0F7gokrceKKopLkOAggggAACCCCQHAIRA15d4mvfvn1y7bXXmnG7OjZ248aN5l7yWVlZcuWVV0q5cuWKTMO58USkC3LjiUhC7EcAAQQQQAABBFJDIOKQBp0INmPGDLn00kvNjSaCWTTYfeGFF4I3uf6cG0+4TswFEEAAAQQQQACBpBKIGPBqbXXlgyVLlsiCBQvMBDa9y9oJJ5xgbitc1Bpu33hC7yz31VdfRayWjlv+66+/Ih7HAQgggAACCCCAAAKJFYgq4NUipqeny9FHH23G7GZmZprXiSp68I0ndNUGXSLt8OHD0rhx40LfeEKHbCxbtixi1fTuc3rjCxICCCCAAAIIIICAtwWiCnhXr14td9xxh4wfP97U5scff5S3335bdM3dIUOGSJkyZYq0lroO78iRI2Xs2LGizw8ePGiurxPs6tevLzruePDgwVKQyXSnnXaa6E+kpAZFtTJFpLKwHwEEEEAAAQQQQCBvgYgBr/ZiduzY0QSVI0aMkAceeMDkdtZZZ8n1118vO3bskBdffDHvK8R5D+vwxhmU7BBAAAEEEEAAgSQXiBjw6g0ntBd10aJF5ha+gwYNMiRXXHGFaI9qz549xbIsSUtLKxIq1uEtEmYuggACCCCAAAIIJI1AxGXJFi9eLE2bNjXBbs5at2jRQnTy1vLly3Pucu016/C6RkvGCCCAAAIIIIBAUgpEDHgbNGhgliXbtGlTLoB3333XjOOtU6dOrn1ubWAdXrdkyRcBBBBAAAEEEEhOgYhDGvSWwnXr1pWLLrpI7rrrLrMagvb6fvrpp/L888+LrotbsmTJItNhHd4io+ZCCCCAAAIIIIBAUghEDHh1pYMPP/xQevXqZYJbrXX37t1N5Tt37izPPPNMkUK4vQ5vkVaGiyGAAAIIIIAAAgi4LhAx4N22bZu54YTebW3+/Pmivbvao6vjeps0aeJ6AcNdIHgd3pUrV5oxxNnZ2VK7du1Cr8Mb7npsQwABBBBAAAEEEPCvQMSA97PPPjM9uqtWrRKdpKY/XkkaeGvwqz+adu/eLVu2bJEqVap4pYiUAwEEEEAAAQQQQCDBAhEnrVWtWtUUcfv27QkuauTLf/DBB9K8efPIB3IEAggggAACCCCAQMoIROzh1TGzOlHs7LPPlssuu0waNmwoGRkZIUA6ma2o0hNPPJHnMmh//PGH6eW99dZbTXFOPvlkue6664qqaFwHAQQQQAABBBBAwIMCEQPe2bNny0cffWSK/s4774StQlEGvL/88ovocmiVKlWSnMuh6XhjHcs7ffp0U86cgXnYwrMRAQQQQAABBBBAIKkFIga87du3N72mXlEYM2aMnHjiiTJs2DDp0aOH9OvXT4oV+7+RGW+++abceeed8uuvv3qluJQDAQQQQAABBBBAIMECYcfwHjx4UPbt25fgooW/vAa39913n0ybNk1eeeUVadu2raxYsSL8wWxFAAEEEEAAAQQQSHmBsAGvrq2bmZkZgrN69Wr57bffQrYl8sUpp5wis2bNMsujaY/vG2+8kcjicG0EEEAAAQQQQAABjwpEHNLglHvEiBFmbOzMmTOdTQl/LFOmjDz33HNy8cUXm8lppUuXTniZKAACCCCAAAIIIICAtwTC9vB6q4iRS6MBr47bbdmypeiqEiQEEEAAAQQQQAABBByBqHt4nRO8+li9enXJaxUJr5aZciGAAAIIIIAAAgi4L5AUPbzuM3EFBBBAAAEEEEAAAb8KEPD6teUoNwIIIIAAAggggEBUAnkOadi0aZMcf/zxgUzWr19v1uMN3ubsnD9/vvOURwQQQAABBBBAAAEEPCUQNuBt1KiRXHLJJSEFPfroo0Ne8wIBBBBAAAEEEEAAAT8IhA14O3XqJPpDQgABBBAoWoG0tDQpXry4axfV9cubN2/uWv5kjAACCHhRIGzA68WCUiYEEEAgFQQsy5KpU6eKBr7xTrqeuq6lTsAbb1nyQwABrwsQ8Hq9hSgfAgiknIAGu3obdRICCCCAQHwE+IsaH0dyQQABBBBAAAEEEPCoAD28Hm0YihUqsHv3btm7d2/oxji9Sk9PlwoVKsQpN7JBAAEEEEAAAa8JEPB6rUUoTy6BHTt2mIC0YsWKufbFY8O2bdtEl9Zr2rRpPLIjDwQQQAABBBDwmAABr8cahOLkFtDe3apVq8r48eNz74zDlvvvv1+WL19OwBsHS7JAAAEEEEDAiwIEvF5sFcpUpAKbN2+WrKwsqVWrVtyvu3r1amnRokXc8yVDBBBAAAEEEIhegIA3eiuOTFKB7du3y4ABA6Rhw4Zxr+FHH30kq1atinu+ZIgAAggggAAC0QsQ8EZvxZH5CKxbt07at28vNWrUyOeogu3atWuXZGdnF+zkKM8qVaqUHHXUUVEeHf1hOiGOhAACCCCAAAKJFSDgTax/0lz9lVdekQYNGkjLli3jXqetW7fKvHnz4p4vGSKAAAIIIIBAaggQ8KZGOxdJLcuXL+9KwLt48eIiKT8XQQABBBBAAIHkFODGE8nZrtQKAQQQQAABBBBA4P8FCHh5KyCAAAIIIIAAAggktQABb1I3L5VDAAEEEEAAAQQQIODlPYAAAggggAACCCCQ1AIEvEndvFQOAQQQQAABBBBAgICX9wACCCCAAAIIIIBAUgsQ8CZ181I5BBBAAAEEEEAAAQJe3gMIIIAAAggggAACSS1AwJvUzUvlEEAAAQQQQAABBHwf8FqWJZs2bZItW7bQmggggAACCCCAAAII5BLwZcC7Zs0aGTBggNSvX19Kliwp1atXlypVqkiFChWkWbNm0q9fP9m1a1euyrIBAQQQQAABBBBAIPUESvityitWrJDWrVtLWlqaZGVlSWZmplSuXNm81l7eZcuWybhx42T8+PEyZcoUadiwod+qSHkRQAABBBBAAAEE4ijgu4B3+PDhpmd38uTJkpGREZZi2LBh0qFDBxk9erQMHjw47DFsRAABBBBAAAEEEEgNAd8NaZgzZ4706NEjz2BXmy09PV169eolkyZNSo1WpJYIIIAAAggggAACeQr4LuBt1aqVzJgxI88KOTumTp0qderUcV7yiAACCCCAAAIIIJCiAr4b0tC1a1fRoHfDhg3SrVs3M0ZXJ6wVK1bMrNSwfPlyGTNmjEycOFF02AMJAQQQQAABBBBAILUFfBfwnnTSSTJv3jzp06eP9OzZUw4fPpyrBdu1aydffvmltGnTJtc+NiCAAAKpKrB161YZOHCgvP7663En2Llzp/m7fMstt8Q9bzJEAAEECivgu4BXK9yoUSOzAsOBAwdk5cqVor262dnZUrt2balbt65ZoqywMJyPAAIIJJuALul42WWXSZMmTeJetd27d8uLL74oBLxxpyVDBBCIg4AvA16n3roGry47puvvFi9e3CxP5uzjEQEEEEAgVECXcyxbtqzoN2XxTgsXLpQSJXz9X0q8ScgPAQQ8JOC7SWtqx40nPPQOoigIIIAAAggggIDHBXz3cZwbT3j8HUXxEEAAAQQQQAABjwn4LuDlxhMeewdRHAQQQAABBBBAwOMCvhvSwI0nPP6OongIIIAAAggggIDHBHwX8HLjCY+9gygOAggggAACCCDgcQHfDWngxhMef0dRPAQQQAABBBBAwGMCvgt43b7xxJ9//inffvttxGbauHGj6CLuJAQQQAABBBBAAAFvC/gu4FVON288sX79+qgCXr2r0J49e7zdupQOAQQQQAABBBBAQHwZ8Drt5saNJ1q3bi36Eym1bNlS6tSpE+kw9iOAAAIIIIAAAggkWMB3k9bUixtPJPhdw+URQAABBBBAAAEfCfiuh5cbTxT83aXjk7/55puCZ5DPmb/88ouUK1cunyPYhQACCCCAAAIIJEbAdwEvN54o+BvlmGOOkVNPPVVq1KhR8EzyOHPChAnSvXv3PPayGQEEEEAAAQQQSJyA7wJevfFEz549JSMjI0+19PR06dWrl4waNUoGDx6c53GptkPHHPfr109q1aoV96qvXbs27nmSIQIIIIAAAgggEA8B343h5cYT8Wh28kAAAQQQQAABBFJHwHc9vNx4InXenNQUAQQQQAABBBCIh4DvAl63bzwRD1TyQAABBBBAAAEEEPCOgO8CXqVzbjyxf/9+WbVqlSxfvlyys7Oldu3aUrduXalSpYrs3bvX/JQuXdo72pQEAQQQQAABBBBAoMgFfDeGV4VGjx4tDRo0kPLly8s111wjGtR26NBBmjVrZoJdPUaHPujkNhICCCCAAAIIIIBAagv4LuD96quvTCBbr149ufvuu2XTpk1y9tlny7PPPpvaLUntEUAAAQQQQAABBMIK+G5IwwsvvCAXXnihTJo0yVTokUcekUGDBknfvn3NjQ/o1Q3bzmxEAAEEEEAAAQRSVsB3Aa/eaS04qE1LS5MhQ4bIoUOHpHfv3qJrzbZr1y5lG5SKI4AAAggggAACCIQK+G5Ig05Mmzp1amgt7FdDhw6Vq6++Wrp06SLz5s3LtZ8NCCCAAAIIIIAAAqkp4LuAVyejffrpp6Y3d/bs2SGt9uqrr8o555wjbdq0kV9//TVkHy8QQAABBBBAAAEEUlPAdwHvVVddJQ899JC8+eab5ie42UqUKCFjx46Vzp07y5IlS4J38RwBBBBAAAEEEEAgRQV8N4ZX22ngwIFmhYZt27blaraMjAx57bXX5KabbpINGzbk2s8GBBBAAAEEEEAAgdQS8GXAq01UqlQpqVmzZp6tdfrpp+e5jx0IIIAAAvEX2LJlizz66KPxz9jO8bjjjpPLL7/clbzJFAEEkl/AtwFv8jcNNUQAAQT8I6B3vNSfRYsWuVJo/WZv/vz50rRpU1fyJ1MEEEhuAQLe5G5faocAAggUiYBlWVK/fn254YYbXLned999J3oNEgIIIFAQAd9NWitIJTkHAQQQQAABBBBAIHUFCHhTt+2pOQIIIIAAAgggkBICBLwp0cxUEgEEEEAAAQQQSF0BAt7UbXtqjgACCCCAAAIIpIQAAW9KNDOVRAABBBBAAAEEUleAgDd1256aI4AAAggggAACKSFAwJsSzUwlEUAAAQQQQACB1BUg4E3dtqfmCCCAAAIIIIBASggQ8KZEM1NJBBBAAAEEEEAgdQUIeFO37ak5AggggAACCCCQEgIEvCnRzFQSAQQQQAABBBBIXQEC3tRte2qOAAIIIIAAAgikhAABb0o0M5VEAAEEEEAAAQRSV6BE6ladmiOAAAII+EXg4MGD8t1338natWtdKfJ5550nxYsXdyVvMkUAgcQLEPAmvg0oAQIIIIBABIENGzbIkCFD5Mgjj4xwZOy7f/jhBxk6dKjcf//9sZ/MGQgg4AsBAl5fNBOFRAABBFJboGTJkvLoo49KZmZm3CHee+890YCahAACySvAGN7kbVtqhgACCCCAAAIIIGALEPDyNkAAAQQQQAABBBBIagGGNHiseSdMmCDr1693pVQ66YOEAAIIIIAAAgikmgABr4dafMaMGXLppZeaHzeKtW3bNjeyJU8EEEAAAQQQQMDTAgS8Hmqeffv2SYsWLaR///6ulGrq1Kmu5EumCCCAAAIIIICAlwUYw+vl1qFsCCCAAAIIIIAAAoUWIOAtNCEZIIAAAggggAACCHhZgIDXy61D2RBAAAEEEEAAAQQKLUDAW2hCMkAAAQQQQAABBBDwsgCT1rzcOpQNAQQQQKBIBDZu3Cg//vijK9eqV6+e1KxZ05W8yRQBBKITIOCNzomjEEAAAQSSVGDp0qXy+eefy/z58+New/3798vixYvFsqy4502GCCAQvQABb/RWHIkAAgggkIQCO3bskMsvv1xuv/32uNduz549csUVV8Q9XzJEAIHYBBjDG5sXRyOAAAIIIIAAAgj4TICA12cNRnERQAABBBBAAAEEYhNgSENsXhyNAAIIIIBATAKHDx+Wr776KqZzoj24bNmy0rJly2gP5zgEUlaAgDdlm56KI4AAAgi4LaCT1vS28QMHDnTlUrqyhN42vm3btq7kT6YIJIsAAW+ytCT1QAABBBDwnMChQ4ckPT1dHn/8cVfKNmTIENm0aZMreZMpAskkQMCbTK1JXRBAAAEEUkpA1w/u2rWr3HPPPXGvtwbSDz/8sNx9991xz5sMEShqAQLeohbneggggAACCMRJYOvWrSbYPf744+OU49/ZLF++XD777DMC3r9JeOZjAQJeHzceRUcAAQQQSG2BtLQ0KVmypNStWzfuEAsWLJDp06dLqVKl4p63Zqjjm3UN5HLlyrmSP5kiECzg+4BX716zefNmKV68uFSuXDm4bjxHAAEEEEAAgQIK6JCGpk2byogRIwqYQ/6nde/eXbZt20bAmz8Te+Mk4MuAd82aNTJy5EgZO3as6PODBw8ajvLly0v9+vWlXbt2MnjwYNHlWkgIIIAAAgggUDAB7UHOyMgo2MkRztK8SQgUlYDvAt4VK1ZI69atRX9RsrKyJDMz0/Ts6ustW7bIsmXLZNy4cTJ+/HiZMmWKNGzYsKgsuQ4CCCCAAAIIRCmgnVWvvfaaVKpUKcozoj9MY4Jbb701+hM4MukFfBfwDh8+XOrbvbiTJ0/O81PnsGHDpEOHDjJ69GjT0xtLK2pA/csvv0Q85a+//pLt27dHPC7WA3bu3GnGTMV6XjTHZ2dny/z582XRokXRHB7TMbrOpHroeK94J/0gc+DAAVfy1rLu2rXLLOvjRtl3795t8ncjby27Lnn0888/myE9+jqeae3atebrRjfKvnr1atm7d69rbap5r1y50nwIjqeJ5rVnzx6zrqobLk5Zv/nmG/Oh3nkdr0f9ilq/Qnaj7OvWrTM2buSt9de/Afq3a9WqVfHiCOSj7xe3/n7p+0X/9rrlou3pVtl1Qpz+DXOr7Pr/qHZOVatWLdAW8XqiHV59+/aVGjVqxCvLQD469ljz7dWrV2BbPJ/88ccf0rhx43hmGZKXDiVxY8x3yEU8+CLNHgNrebBceRZJe3d79uwpvXv3zvMY3fHWW2/JqFGjRBfljiV9/fXX8uyzz0Y8RQNHPe7888+PeGy0B+gfluuvv94EAtGeE8txGkzrV1M6wSHeSYNG/alZs2a8szZBnf6n4dbEBp00ob8GFSpUiHvZ9T9S/Y/ajby1sPofnQ7d0THs8U663JG6uPEfhgYA+p+GW8OO9L2uPTxu5K9tqj9uzRnQAKZixYrxbk6TnwYY2qZVq1aNe/7apvrB163fU/37UqJECVcmUOnMUd0HAAAYpklEQVTfXv0b40bgpXdZ0/ejW38D9O+X9pS68X50bprhVtnd/D9JP7Brm7rxf5L+Dv3www/SoEGDuP8eaYa//fabHHXUUXLEEUfEPX99r2jHYceOHeOet9cz9F3Aq2sN6n/Er7/+er621113nenJ+OCDD/I9jp0IIIAAAggggAACyS3guyENusB2q1atZMOGDdKtWzczRrdKlSpSrFgx8/Wlrhs4ZswYmThxohn2kNzNR+0QQAABBBBAAAEEIgn4rodXK/Tnn39Knz59ZNq0aaJfF+VMukrD/fffz73Fc8LwGgEEEEAAAQQQSEEBXwa8Tjvp2EidmKK9ujp+rHbt2mYgtvb4khBAAAEEEEAAAQQQUAFfB7w0IQIIIIAAAggggAACkQSKRTqA/QgggAACCCCAAAII+FmAgNfPrUfZEUAAAQQQQAABBCIKEPBGJOIABBBAAAEEEEAAAT8LEPD6ufUoOwIIIIAAAggggEBEAd+twxuxRhyAQIwCDz/8sLlTmRt3KoqxKBweJ4E5c+ZIZmamlC9fPk45kk2iBWbOnClNmjRx5e5Tia5bql5f71b25JNPSvXq1VOVgHoXoQABbxFicylvCowbN07q169vfrxZQkoVq8CECRPMDWpq1aoV66kc71GBTz75xNy6mGUnPdpABSjWhx9+KPfeey8BbwHsOCV2AQLe2M04I8kE6tWrJzfffLNcdNFFSVaz1K3O7NmzpV+/fiboTV2F5Kq53mjowQcflKZNmyZXxVK4Nh9//LFUqFAhhQWoelEKMIa3KLW5FgIIIIAAAggggECRCxDwFjk5F0QAAQQQQAABBBAoSgEC3qLU5loIIIAAAggggAACRS5AwFvk5FwQAQQQQAABBBBAoCgFCHiLUptrIYAAAggggAACCBS5AAFvkZNzQQQQQAABBBBAAIGiFCDgLUptroUAAggggAACCCBQ5AIEvEVOzgURQAABBBBAAAEEilIgzbJTUV6QayHgNYHly5eL3r2pXLlyXisa5SmgwJIlS0TvslamTJkC5sBpXhNYtGiR6E1iMjIyvFY0ylNAgd9//10aNmwo6enpBcyB0xCIXoCAN3orjkQAAQQQQAABBBDwoQBDGnzYaBQZAQQQQAABBBBAIHoBAt7orTgSAQQQQAABBBBAwIcCBLw+bDSKjAACCCCAAAIIIBC9AAFv9FYciQACCCCAAAIIIOBDAQJeHzYaRUYAAQQQQAABBBCIXoCAN3orjkQAAQQQQAABBBDwoQABrw8bjSIjgAACCCCAAAIIRC9AwBu9FUcigAACCCCAAAII+FCAgNeHjUaREUAAAQQQQAABBKIXIOCN3oojEUAAAQQQQAABBHwoQMDrw0ajyIUTsCyrcBlwticECtqOBT3PE5VO8kLQNsnZwAcOHJBDhw7FVDneCzFxcXAUAgS8USBxSHIIzJ49W7p16yaVKlWSzMxMeeSRRyJWbPz48XLsscfm+tHtpMQIFKQdd+3aJQMGDJCjjz5aKleuLJdffrn89ddfiakAV80lUJA25XczF6MnNyxbtkxq164tkyZNilg+fk8jEnFAIQRKFOJcTkXANwJ79uyRLl26yOmnny5ff/21zJ07V2655RYpVqyYPPDAA3nW47vvvhM9t1evXiHHNGjQIOQ1L4pGoKDteP/998vEiRPl+eefl5IlS8ptt90m7dq1k1mzZklaWlrRFJ6rhBUoaJvyuxmW01MbNdjt2LFj1B8u+T31VPMlX2Hsrw1ICCS9wKBBg6zy5ctbe/fuDdR18ODBVtWqVa19+/YFtuV8ct5551l2sJtzM68TJFCQdpw3b55lf7CxPvroo0CpFy5cqONaLLvXKbCNJ4kRKEibakn53UxMe0V71eeee8464ogjrMaNG5vftQkTJuR7Kr+n+fKwMw4CDGlIvs8w1CiMwBdffCEdOnSQUqVKBfZ26tRJNm/eLD///HNgW84n9h9had68udl8+PDhnLt5XcQCBWnHr776yvTqavs7qUmTJnLMMcfIZ5995mziMUECBWlTLSq/mwlqsCgvO3ToUOnbt6/5ZiWaU/g9jUaJYwojQMBbGD3O9Y3An3/+KXXq1Akpr/N6/fr1IdudF7p906ZN8vvvv8u5555rguWTTz5Zvv32W+cQHotYoCDtqOdUq1bNBL3BxdVxhRs2bAjexPMECBSkTfndTEBDxXhJHZf92GOP5fq9yysbfk/zkmF7vAQIeOMlST6eFtixY4dUqVIlpIwVK1Y0r/MKenScr6YZM2bIVVddJQ899JAJgDX41d4lUtELFKQd9RydqJYz6eTFvNo+57G8dk+gIG3K76Z77RGvnPVDZiyJ39NYtDi2IAJMWiuIGud4VmDbtm2iP05KT083Pbv6mHNykvM6OzvbOTzkUb/2fvHFFyUrK0uc4Lh3795mxvGQIUNk3LhxIcfzwn2BgrSjnqOTE3MmbX9dLomUWIGCtCm/m4ltMzeuzu+pG6rkGSyQ+3+B4L08R8BnAqNGjRJdQcH5ad++valBzZo1ZevWrSG1cQLjcuXKhWx3Xhx11FGiAa4T7Or2GjVqyNlnn21WeXCO47HoBArSjuHO0RLr+8GeyFh0hedKYQXCtQ+/m2GpknpjuPeBVpjf06Ru9iKtHD28RcrNxdwW0CVw6tatG7iMfm2tSf+Y5hyru27dOrOvYcOG5jHnPytXrpS1a9fKGWecEbJLJ75pbwSp6AUK0o61atUyQ1F04fvixYsHCq3vh3POOSfwmieJEShIm/K7mZi2cvOq/J66qUveKkAPL++DpBJo1qyZXHvttYGfzp07m/rpmqu68Hnw3X50hr728J1yyilhDV599VVp1aqVmbTmHKA9T1OmTJFTTz3V2cRjEQoUpB3t5atk9+7dMm3atEBJly5dKr/99puZjBjYyJOECBSkTfndTEhTuXpRfk9d5SVzFYjD0mZkgYDnBeyeWsvulbXsm01YW7Zssezgx7InMllPP/10oOxvvvmmdeWVV1r2Qvhm2x9//GGVLl3aOv/88605c+ZY9gQnq3v37lZGRoa1YMGCwHk8KTqBaNpx8eLFph2nTp0aKFjbtm2t448/3tI2Xb16tWVPPLRat25t2UvNBY7hSWIECtKm/G4mpq0KctVVq1aFXYeX39OCaHJOYQQIeAujx7m+EtCFz/VGE/o5T4PdG2+80Tp48GCgDv379zf77NnCgW16jj2W12zX8xo1amTZd3gK7OdJ0QtEascff/zRtNfLL78cKJwGuXZvvdleokQJ8yFGbz5B8oZAQdqU301vtF2kUuQV8PJ7GkmO/fEWSNMMtaeXhEAqCOjbfcWKFWacrx34RFVlPccOmMx6kjppjZR4gYK0o5ZabzSiKzaEW6Ys8bVK7RIUpE353UzO9wy/p8nZromuFQFvoluA6yOAAAIIIIAAAgi4KsCkNVd5yRwBBBBAAAEEEEAg0QIEvIluAa6PAAIIIIAAAggg4KoAAa+rvGSOAAIIIIAAAgggkGgBAt5EtwDXRwABBBBAAAEEEHBVgIDXVV4yRwABBBBAAAEEEEi0AAFvoluA6yOAAAIIIIAAAgi4KkDA6yovmSOAAAIIIIAAAggkWoCAN9EtwPURQAABBBBAAAEEXBUg4HWVl8wRQAABBBBAAAEEEi1AwJvoFuD6CCCAAAIIIIAAAq4KEPC6ykvmCCCAAAIIIIAAAokWIOBNdAtwfQQQQAABBBBAAAFXBQh4XeUlcwQQQAABBBBAAIFECxDwJroFuD4CCCCAAAIIIICAqwIEvK7ykjkCCCCAAAIIIIBAogUIeBPdAlwfAQQQQAABBBBAwFUBAl5XeckcAQQQQAABBBBAINECBLyJbgGujwACCCCAAAIIIOCqAAGvq7xkjgACCCCAAAIIIJBoAQLeRLcA10cAAQQQQAABBBBwVYCA11VeMkcAAQQQQAABBBBItAABb6JbgOsjgAACCCCAAAIIuCpAwOsqL5kjgAACCCCAAAIIJFqAgDfRLcD1EUAAAQQQQAABBFwVIOB1lZfMEUAgmQUOHDgghw4dyrOKlmXluS+vHQcPHpTly5eL5p3otGrVKlm/fn2ii8H1EUAAgUILEPAWmpAMEEAgkQK//fabpKWlydtvv12kxVi2bJnUrl1bJk2alOu6b7zxhpxzzjlSpkwZOe2002TatGkhx9SoUcOUWcs9d+5cs2/p0qVy1llnyRFHHCENGjSQsmXLykknnSQ//vhjyLmVKlWSf//73yHb4v1CA9369evLUUcdJc2bN49L9ocPH5b/+Z//kX379hUqv06dOgXsXnjhhULlxckIIJA6AgS8qdPW1BQBBOIkoMFux44d5a+//sqV4zfffCM33nijdOnSRb7//ns55ZRTpEOHDjJv3ryQY2+55RaZP3++NG7c2ORz6qmnyp49e+TJJ5+UH374QUaNGiXlypWTNm3ayHfffRc49+qrr5bjjz8+8NqNJy+++KKsXbvWBOo//fRTXC7x3nvvyc0335xvj3g0F1IXdSMhgAACsQiUiOVgjkUAAQRSXUB7Ke+++26pW7duWIqbbrpJsrKy5NZbbzX79fj//ve/8swzz8grr7wSOEd7eZs2bWpejxs3TrZt2yZvvfWWNGnSxGw7/fTT5fLLL5d69eqZ7a1atQrkF8jEpScbNmwwZdNgO15Je3jjkbTXmYQAAgjEKkAPb6xiHI8AAr4U+Pbbb+Xss8+WChUqyNFHHy0DBgyQ/fv3h9RlwoQJctFFF0nlypWlXbt2Mn36dNGeV+3RddLQoUOlb9++MnHiRGdT4HH16tWiQywuu+yywDZ9ol/DhzveOUjH7ep43927dzubzGO1atXkk08+Ee3VdZIOlXj99dfNy9tvv92UT8uY82fHjh3mmF27dpmeVR2iULVqVencubOsWLHCyS7XowbZH374oSxatMjkOX78eHOMDrnQfZqP9jxrz/WYMWNCzv/555/lqquukurVq0vbtm3lP//5j2igq3V48MEHzbGtW7eWN9980zzX8c8jRowwQb7mqcM/nOvpAfpBQQP9L7/80nzA0AB8y5Yt5lz+QQABBGIRIOCNRYtjEUDAlwK//PKLCcB0/Ovo0aNNwKrjPzWAc9KsWbPMMATtedVjTjjhBLn44otFzw0edzp79mx57LHHpGTJks6pgcfFixeb53Xq1Als0yf6etOmTSb4C9nx/y/at29vglEd+qCBoY7bdXpENfDWQN1JOuZXe2A1XXDBBXLdddeZn+uvv1569OhhAu5ixYqZ8cMaRJ933nkmiOzZs6e89NJLsnXrVmnRokWegaMepz3PNWvWFA2o1WH79u0mwNU6DB482ASppUqVku7duweGF+iHAv2woMMyXn75ZTOG+a677jJjq4877ji58MILTZl1KIcGy5oGDRpkPnhokPzOO+/IySefbNrgtddeM/v1ujqkQoeIqIF+ENEfEgIIIBCzgP0HkYQAAgj4VmDhwoW6FIJl9zbmWQd7Mphl94BadhAZOMYO/sx5X331ldlmB1SWHegF9usTO4Azx+g1ciZ7YpfZZ/cKB3Z9/PHHZpvdOxrYpk/sYM5s37hxo9lu94BaQ4YMCTnG7hm27Ali5jitjx2cW127drXsADfkuIoVK1r/+te/QrY5L6688kpL89ayaXKua/cuO4dYdkBq2RPjrPvuuy+wLeeTG264wbJ7WwObP//8c+uYY46x7J7hwDY7uDdltcf7mm3/+Mc/rGOPPdaye20Dx9hjdi111aTto/Wye5zN65UrV1rp6em56nL++eebOti975ba6jl2r7o5J/gf3f78888Hb+I5AgggkKcAPbz2X00SAggkr4D9109mzpwpdvBoZvc7NdVeVU06IUyP0Z5b7TENTk6vZPC2/J6XKPF/0yJ09YVwKb+lxuxgUbSXWSdkDR8+3PSCvv/+++Yr/U8//TRcdiHbdKjFBx98IHqOM75Yh3HoEA4dLqD11B+tZ7NmzUImwoVkFOaFWv3xxx9m1YadO3eaXlc7CBbtSdYeXU1z5swxQzl0m5N0SIMOCwmX9Pjs7Gzp1q1byG7t5bY/GMiSJUsC2+M5ljiQKU8QQCClBJi0llLNTWURSD0B/fp/7969Zgmx4NprUJiZmWlWI9AVCTSQO+OMM4IPkTPPPDPkdaQXtWrVMofoBLTgpMMINJUvXz54c9jnOpxAf/r37y9//vmnCSJ1aMGll14a9njdaPcsy8CBA2XkyJEhwx90PV8dFqDLneVMOhY3lqST7nSFBB3LW7p0aROI6wcF/dHhFzrmV4dBRJt0CIR+MHDMnPOc4FbbxElOAO+85hEBBBCIVeDvj+KxnsnxCCCAgA8EdMyn9jqGm+ykE7oaNmxoxs9q8OUEpk61NFiMJTkB37p160JO05s36IQx7WkNlzQgveaaa3LtatSokdx5551m0lxeN4BYsGCBOVfH3jorQzgZValSxazyoJPDNCgN/tHANdqkq0toOXTC24wZM8yKEhpka9KAV301mM8Z6OskvJwWzjV1Qp6em9NcP3ho0nZxUvHixZ2nPCKAAAIFEiDgLRAbJyGAgF8EdHKZrsowefLkkCLr5C/96vzEE0+UjIwMM4Qg59ABXV0glqS9ldo7+9lnn4Wcpqs/nHvuuSHbgl/o6gR67XBBqA5x0EBaJ9PlTBrE6woQOhxClz/LmbQsuiKDToLTgF5/NMjU4R1PPPFEzsPzfK0rTOjEMz2nZcuWYo+9NUMYNC/nTnO6X1dVCE5PPfWUmYimk/702po06NakZdOUs130tQbP9WPsgTaZ8Q8CCCCQhwBDGvKAYTMCCPhLYOrUqbmW9dIa9O7dW+6//37p1auXaACmqxnYE6ZEe0R1VQBdJkvTQw89ZIJHHbfbp08fMxZW18eNNWkvq/aG2pOvzI+OY7UnvZn88spLy6h3itPAV8/VVRR0GIaOk9XVJB555JFAwBichz1Jzdz6V+ulqxloAOokDSi1HjoeWFdGePjhh80yY08//bQpyz333OMcGvFRV0jQHt2vv/7aeOmYaGfsrTOGV/PT4Fvzv/baa03w+9xzz5mhGbqig9O7rWsNq7F+0LjkkkukX79+JqDXOusHA12mTA1ICCCAQFwF7D+QJAQQQMC3As4qDfYfRjOjP/jR7lU09bJ7FS07ELPsoMscYwdglh10Wc6qCU7l7buBWXYQbNk9jJa9nJf17LPPmuPtsbDOIYHHcKs06E57TV3LXo7LrECgZbFvJGHZtxoOnKdPwq3SoCsg2EuQmWs7dbBvOmHZ43JDznVWabAnwOWqr3OePtrr2Zrz7ElqgdUf1EPrF2l1g5yrNNjDDix7OINlj9217Il5lt3jbNmBqzG0J5kFymcHuCHl13OcVRk0D/tWyabMunqDps2bN1v2kmSWPWTBbFeX4NUjnFUanFUnAheyn2gdI9Uj+HieI4BAaguY/w3sPxwkBBBAIOkF7D/35it+HXqgwxiCk940QSex6bhXJ+kNErRnWMfyli1b1tkc1aN+ja9DJsLdGUyHJ2hPsE40C5e0B7pMmTJm3G+4/QXZpuNrtUzOOOOC5KG9uTqMIr9JZGpsf0AQXfPYDs5zXUbLoZbOihZ6gJZLx/o2aNAg1/F5bdAhEnbAa3qx8zqG7QgggIAjwJAGR4JHBBBIegENkurnMTZUg08NUPVre13KS1d3ePTRR82qB7EGuwqpX+OHC3ajQS7oefnlHS74zO/4cPs0CNef/JIa5xe4hiuHWuV3Tn7XYx8CCCAQjQCT1qJR4hgEEEh6AftmDuZWw7qagq4QoL3A9g0a8h17WxgUHVOrvZw6eY4UvYDetjm4dzj6MzkSAQRSWYAhDanc+tQdAQRyCdjjRQM3Z7DH0ObaH48N9l3VxB7ra7LSFSS0h5MUnYCuOrFjxw5zsN6ymVsNR+fGUQikugABb6q/A6g/AggggAACCCCQ5AIMaUjyBqZ6CCCAAAIIIIBAqgsQ8Kb6O4D6I4AAAggggAACSS5AwJvkDUz1EEAAAQQQQACBVBcg4E31dwD1RwABBBBAAAEEklyAgDfJG5jqIYAAAggggAACqS5AwJvq7wDqjwACCCCAAAIIJLkAAW+SNzDVQwABBBBAAAEEUl2AgDfV3wHUHwEEEEAAAQQQSHIBAt4kb2CqhwACCCCAAAIIpLoAAW+qvwOoPwIIIIAAAgggkOQCBLxJ3sBUDwEEEEAAAQQQSHUBAt5UfwdQfwQQQAABBBBAIMkFCHiTvIGpHgIIIIAAAgggkOoCBLyp/g6g/ggggAACCCCAQJILEPAmeQNTPQQQQAABBBBAINUFCHhT/R1A/RFAAAEEEEAAgSQXIOBN8gameggggAACCCCAQKoLEPCm+juA+iOAAAIIIIAAAkkuQMCb5A1M9RBAAAEEEEAAgVQX+F9LQ9JAAT81PgAAAABJRU5ErkJggg==" />

<!-- rnb-plot-end -->

<!-- rnb-chunk-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuc2V0LnNlZWQoMTAwKVxuIyBjbHVzdC56ZWlzZWwgPC0gcXVpY2tDbHVzdGVyKHNjZS56ZWlzZWwpIFxuIyBzY2UuemVpc2VsIDwtIGNvbXB1dGVTdW1GYWN0b3JzKHNjZS56ZWlzZWwsIGNsdXN0ZXI9Y2x1c3QuemVpc2VsLCBtaW4ubWVhbj0wLjEpXG5saXZlcl9zY2UgPC0gbG9nTm9ybUNvdW50cyhsaXZlcl9zY2UpXG5hc3NheU5hbWVzKGxpdmVyX3NjZSlcbmBgYFxuYGBgIn0= -->

```r
```r
set.seed(100)
# clust.zeisel <- quickCluster(sce.zeisel) 
# sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clust.zeisel, min.mean=0.1)
liver_sce <- logNormCounts(liver_sce)
assayNames(liver_sce)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiWzFdIFxcY291bnRzXFwgICAgXFxsb2djb3VudHNcXFxuIn0= -->

[1]  

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


### Feature selection


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubGlicmFyeShzY3JhbilcbmRlY19saXZlciA8LSBtb2RlbEdlbmVWYXIobGl2ZXJfc2NlKVxuXG4jIFZpc3VhbGl6aW5nIHRoZSBmaXQ6XG5maXRfbGl2ZXIgPC0gbWV0YWRhdGEoZGVjX2xpdmVyKVxucGxvdChmaXRfbGl2ZXIkbWVhbiwgZml0X2xpdmVyJHZhciwgeGxhYj1cIk1lYW4gb2YgbG9nLWV4cHJlc3Npb25cIixcbiAgICB5bGFiPVwiVmFyaWFuY2Ugb2YgbG9nLWV4cHJlc3Npb25cIilcblxuaHZncyA8LSBnZXRUb3BIVkdzKGRlY19saXZlciwgbj0zMDAwKVxuYGBgIn0= -->

```r
library(scran)
dec_liver <- modelGeneVar(liver_sce)

# Visualizing the fit:
fit_liver <- metadata(dec_liver)
plot(fit_liver$mean, fit_liver$var, xlab="Mean of log-expression",
    ylab="Variance of log-expression")

hvgs <- getTopHVGs(dec_liver, n=3000)

Dim reduction

liver_sce <- scater::runPCA(liver_sce, subset_row=hvgs, ncomponents=30)
reducedDim(liver_sce, "PCA") <- reducedDim(liver_sce, "PCA")[,1:11]
plotPCA(liver_sce, colour_by="Condition", ncomponents=3)

liver_sce <- runUMAP(liver_sce, dimred="PCA", ncomponents=2)

scater::plotUMAP(liver_sce, colour_by="Condition", point_alpha=1,  point_size=0.8) 
scater::plotUMAP(liver_sce, colour_by="Sort", point_alpha=0.3,  point_size=0.5)
scater::plotUMAP(liver_sce, colour_by="Patient", point_alpha=0.3,  point_size=0.5)